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FOREWORD 

This final report was prepared by personnel of the Computa- 
tional Mechanics Section of Lockheed's Huntsville Engineering 
Center, It constitutes final documentation of efforts performed 
under Contract NAS8-36090 for NASA-Marshall Space Flight Center. 

The NASA-MSFC Contracting Officer's Representative for this 
research study was Dr, P.K. McConnaughey, ED32. 
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1 * INTRODUCTION 


The objective of this research effort, as defined by NASA-Marshall 
Space Flight Center was two-fold: (1) to numerically simulate viscous 

subsonic flow in a proposed elliptical two-duct version of the fuel side Hot 
Gas Manifold (HGM) for the Space Shuttle Main Engine (SSME), and (2) to 
provide analytical support for SSME related numerical computational 
experiments, being performed by the Computational Fluid Dynamics staff in 
the Aerophysics Division of the Structures and Dynamics Laboratory at 
NASA-MSFC. Numerical results of HGM were calculations to complement both 
water flow visualization experiments and air flow visualization experiments 
and air experiments in two-duct geometries performed at NASA-MSFC and 
Rocketdyne. In addition, code modification and improvement efforts were to 
strengthen the CFD capabilities of NASA-MSFC for producing reliable 
predictions of flow environments within the SSME. 

Full three-dimensional laminar and turbulent SSME /HGM computations were 
performed on a sparse grid model of the Phase 11+ two-duct version for the 
fuel side manifold employing the Lockheed explicit PAGE Navier-Stokes code. 
The approach, methodology, and results of this effort are documented in an 
earlier interim technical report submitted in March 1986 (Ref. 1). These 
will not be reiterated here, and the reader is referred to that previous 
document for those details. 

Emphasis in this report is placed on the second of the aforementioned 
objectives. The direction of these efforts was to obtain and modify, for 
application to SSME internal flow analyses, state-of-the-art incompressible 
full three-dimensional Navier-Stokes codes. The approach taken is described 
in Section 2. Sample results for both laminar and turbulent computations 
are presented in Section 3. 
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2 . TECHNICAL APPROACH 


2 . 1 GENERAL 


Previous computational fluid dynamics results, reported in the interim 
report for this contract (Ref. 1), were performed with a finite difference 
code employing an explicit solution algorithm. Steady state was obtained 
from a trial initialization by performing successive iterations in time 
until all transients have evolved away. The size of the time step used in 
this procedure is a strong function of the density of points in the grid. 

The higher the density, the smaller the allowable time step. For this 
reason, relatively coarse grids were used, even for regions near the solid 
walls . 

Experience has dictated that much larger nodal densities near the walls 
are desirable for more accurate computational predictions. This precludes 
using an explicit code because of the unnecessarily large number of time 
steps required to obtain a steady state solution. The implicit code INS3D, 
developed at NASA-Ames (Ref. 2), was one of two codes chosen for application 
to additional SSME related computations. The second implicit code consid- 
ered was the FDNS3D incompressible code developed by Y.S. Chen (Refs. 3 and 
4). The implicit solution algorithm incorporated into these codes allows 
for much larger time steps even for grids of larger nodal densities. A 
brief description of the approach and solution methodology for each of these 
codes is now presented for completeness. 

2.1.1 INS3D Code 

The INS3D code solves the three-dimensional incompressible Navier-* 
Stokes equations in primitive variables. An implicit finite difference 
operator is used in a general curvilinear coordinate system. The solution 
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procedure uses the standard approximate factorization scheme. The pressure 
field solution is based on the concept of adding a time-like pressure term 
into the continuity equation via an artificial compressibility factor. This 
approach was first introduced by Chorin (Ref. 5) and later adopted by Steger 
and Kutler (Ref. 6) using an implicit approximate factorization scheme by 
Beam and Warming (Ref. 7). It is from these earlier developments that INS3D 
evolved (Ref. 2). 

Values of the artificial compressibility factor are bounded in order 
not to influence the steady state mass conservation. In the INS3D method- 
ology mass conservation is of crucial importance if a stable solution is to 
result. Since the continuity equation is modified to obtain a hyperbolic- 
type equation, pressure waves of finite speed will be introduced. The speed 
of propagation of these pressure waves depends on the magnitude of the com- 
pressibility parameter. When the pressure waves travel through a given 
location a pressure gradient is created there. Near boundaries, the viscous 
boundary layer must respond to this pressure fluctuation. To accelerate 
convergence and avoid slow fluctuations it is desirable that the time 
required for pressure waves to propagate through the region of interest be 
much less than the time needed for the boundary layer to fully adjust itself. 
This condition provides for a lower bound on the artificial compressibility 
factor. The upper bound on this factor comes not from the physics but from 
the effects of the approximate factorization of the governing equations. 

When the finite difference form of the equation is factored, higher order 
cross-differencing terms are added to the left-hand side of the equation. 
These added terms must be made smaller than the original terms everywhere in 
the computational domain. This condition results in an upper bound on the 
compressibility factor. 

It is well known in the computational fluid dynamics community that the 
approximate factorization schemes which employ alternating direction type 
implicit methods have stability problems in three dimensions (Refs. 8 through 
11). The INS3D code satisfactorily overcomes this difficulty by providing 
second and fourth order smoothing terms to the algorithm to ensure stability 
without adversely affecting mass conservation. 


2-2 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 


LMSC-HEC TR F225988 


2.1.2 FDNS3D Code 


The FDNS3D code solves the steady or unsteady Navier-Stokes equations 
in three-dimensional curvilinear coordinate space for both incompressible 
and compressible flow. In the discretization of the steady transport equa- 
tions, the code employees a combined second and third order upwind differenc- 
ing scheme to approximate the convective terms. Second order central differ- 
encing is used for the remaining terms in the transport equation. For time- 
dependent flow problems, a second order backward differencing scheme is used 
for the temporal term. A SIMPLE type velocity-pressure correction solution 
algorithm is used to couple the continuity and momentum equations. Compress- 
ibility effects are considered in the formulation of the pressure correction 
treatment. Velocity components, pressure, density, and temperature are 
updated using the solution of the correction equation. An alternate direc- 
tion line relaxation matrix solver is applied to solve the system of linear- 
ized algebraic equations. A standard k-e turbulence model and an extended 
version which employs an improved transport equation for the turbulent kin- 
etic energy dissipation rate are included in the FDNS3D code. 

Artificial fourth order smoothing techniques are used in the transport 
equations to obtain strong coupling between the velocity and pressure fields 
in many approaches. To ensure stability of the numerical solution, grid 
staggering between the velocity vectors and pressure nodes are used in a 
usual SIMPLE type algorithm. This method, however, has the drawback that 
the velocity components and the pressure are solved using different control 
volumes. This also requires a lot more bookkeeping in the coding. To ensure 
velocity-pressure coupling for the non-staggered grid system, FDNS3D employs 
an explicit dissipation (smoothing) term which is added to the pressure 
correction equation. This term is limited to a very small contribution to 
the correction equation so that mass conservation is preserved. Both 
approaches, staggered (Ref. 3) and non-staggered (Ref. 4), are coded as 
different versions of FDNS3D. The non-staggered version is used in the 
present analysis. 
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2.2 LAMINAR COMPUTATIONS 

The viscous flows internal to the SSME/HGM involve complicated passages 
which tend to create regions of separation and swirling secondary flow 
patterns. In order to evaluate the incompressible codes on computational 
grids which would possess these characteristics two simple geometries were 
generated. The Lockheed algebraic grid code was employed to develop the 
two-dimensional grid shown in Fig. 2-1 and the three-dimensional grid shown 
in Fig. 2-2. The turnaround duct (TAD) geometry presented in Fig. 2-1 has 
the same dimensions as the Phase 11+ version of the SSME/HGM TAD. Flow 
enters the lower part of the duct from the right, experiences a 180 deg turn 
and exits from the upper part of the duct, flowing toward the right. Adverse 
pressure gradients produced on the downstream side of the turn will tend to 
cause the boundary layer to separate. The three-dimensional flow around the 
turn in the square duct geometry presented in Fig. 2-2 will exhibit a 
swirling secondary flow similar to that exiting the right or left transfer 
ducts of the SSME/HGM. Because of symmetry considerations, no such swirl 
will be produced in the TAD configuration, even in three dimensions. 

Both INS3D and FDNS3D were used to solve identical problems on these 
two geometries. Results of these computations are presented in Section 3.1. 

2 . 3 TURBULENT COMPUTATIONS 

Incorporated into the FDNS3D code for turbulent applications is a k-e 
subroutine with stand and wall function treatment. This was used for the 
two-dimensional case using the TAD geometry. The INS3D code must be supplied 
with an external subroutine for computing the turbulence eddy viscosity 
array. For the two-dimensional TAD geometry, a standard Baldwin-Lomax calcu- 
lation was coded and incorporated into the code. 

For three-dimensional applications, especially to layer computations 
such as primary SSME/HGM geometries, INS3D was chosen to be the primary code. 
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Fig. 2-1 TAD Computational Grid (33 x 117) for Two-Dimensional Test 
Solutions Using Both INS3D and FDNS3D 
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Fig. 2-2 Square Duct Computational Grid (54 x 31 x 27) 
for Three-Dimensional Laminar and Turbulent 
Test Solutions 
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The structure of the finite volume methodology of FDNS3D causes it to be 
about three times more storage intensive as INS3D. On the NASA-MSFC 
Cray-XMP, INS3D can be run as a single zone computation on a geometry 
consisting of no more than 165,000 nodes. FDNS3D can perform the same task 
but on a single zone of one-third as many grid points. Because of this 
limitation and because of the fact that the k-e subroutine in FDNS3D 
seemed to produce consistent and reliable turbulent results (Ref. 12), it 
was decided that this treatment be incorporated into INS3D as a separate 
subroutine. 

The approach for adding the k-e computation to INS3D was to have INS3D 
compute the primative variables, pressure and velocity components, and have 
the turbulence eddy viscosity array updated using two modified subroutines 
from FDNS3D. The subroutines solve the transport equations for turbulent 
kinetic energy, k, and its dissipation rate, e. Thus, the solution proce- 
dures are decoupled. However, since the k-e solution strategy is itera- 
tive at each time stem in INS3D, the decoupling produces a minimal effect on 
the solution. 

In this analysis, the turbulence k-e model after Launder and Spalding 
(Ref. 13) is employed. The modeled kinetic energy and dissipation equations 
are given by 
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where oc is the turbulent kinetic energy dissipation Prandtl number, and 

C and C are model constants. Recommended values for the constants 
1 2 

in the above equations are (Ref. 14): 

C = 1.44 

C =1.92 
2 

C y = 0.09 

o s 1,0 

k 

a ss 1.3 

c 

Much study has been made concerning the values of these constants (Ref. 
15). The values above have produced the most consistent reliable results 
for flows involving streamline curvature and recirculation regions. 


Boundary conditions treatment for the k-e calculation is as follows: 

• Inlet Plane - The strengths of k and e were not known and so were 
determined using 

kin = (turbulence intensity) • (u^ n ) 2 

and 

c in = (k in ) 1 - 5 /(XL), 

where X and L represent a constant and a characteristic length, 
respectively. Sometimes XL is called a turbulent length scale. 

In this analysis, the characteristic length L, length scale X, and 
turbulent intensity are defined AS AY, 0.01 AND 0.0003, respect- 
ively. Parametric study on these variables in a single and multiple 
jet configuration was reported by Bai (Ref. 15). Different flow 
patterns and recirculation zone sizes are obtained near the inlet by 
varying these parameters . The above values were chosen because the 
effect is mainly limited to the inlet regions, as is the effect of 
varying the turbulent model constants C x and C 2 . 

• Exit Plane - For the exit plane, and the symmetric boundary at the 
center of the square duct, extrapolation enforcing normal gradients 
to zero was used. 
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• Solid Walls - Near the solid wall region the flux through the wall 
side interface is set to zero and the diffusive flux term on the left 
side of the equation is transferred into the source term on the right- 
hand side as a false source. This false source term is obtained from 
a "wall function" which is largely based on experimental data (Refs. 

13 and 16). Recently more elaborate wall functions involving the 
division of the near wall region into two or three layers were pro- 
posed. The concept of the wall function is based on the assumption 
that the total stress is equal either to the Reynolds stress for 
the core or to the laminar stress for the viscous sublayer. In the 
current analysis, this wall boundary treatment is applied only to the 
turbulent kinetic energy and its dissipation rate equation. 
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3 . RESULTS 


3.1 LAMINAR COMPUTATIONS 


Laminar results were obtained in the two-dimensional TAD geometry shown 
in Fig. 2-1 using both INS3D and FDNS3D for Reynolds number of 500 and 1000. 
In all cases an essentially converged solution was obtained in 1000 itera- 
tions. Results at a Reynolds number of 1000 are presented in Figs. 3-1 and 
3-2. The two solutions are essentially the same. Both exhibit separation 
on the downstream side of the turn occurring at the same point on the wall 
and both shown the same pressure distribution on the inner and outer walls. 

Since an axisymmetric option is available in FDNS3D, the TAD geometry 
was used to perform such a compilation at the same Reynolds number as the 
two-dimensional case. These results are given in Fig. 3-3. Separation of 
the boundary layer on the inner wall is again observed. It occurs, however, 
much earlier in the turn. This is of course due to the diffusive effect in 
the axisymmetric geometry as the flow traverses the turn. 

A three-dimensional computation was performed with both codes on a 
square duct grid like the one shown in Fig. 2-2 but with the cross-plane 
nodal distribution shown in Fig. 3-4. The difference being that the grid is 
stretched to the solid walls but not to the symmetry plane. Like the two- 
dimensional results these solutions were essentially the same. Thus, only 
the INS3D results are presented, and these are shown in Figs. 3-5 through 
3-7. Again, only 1000 iterations were needed to obtain a converged solu- 
tion. This solution exhibits the characteristic secondary swirl which must 
be produced as the fluid negotiates the 90 deg turn. It also shows no 
separation along the inner wall at the computed Reynolds number of 790. 
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Fig. 3-1 Velocity Vectors, Velocity Magnitude Contours, (Middle) 
and Static Pressure Contours for Re = 1000 in Two- 
Dimensional TAD Using INS3D 
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Fig. 3-2 


Velocity Magnitude Contours, (Bottom) and Static Pressure 
Contours for Re = 1000 in Two-Dimensional TAD Using FDNS3D 
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Fig. 3-4 Square Duct Computational Grid (54 x 41 x 31) for Initial 
Three-Dimensional Laminar Test Solution 
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Fig. 3-7 Velocity Vectors, Velocity Magnitude Contours (Middle) and Static 

Pressure Contours at Square Duct Mid Plane (Left) and Three-Quarter 
Plane for Laminar Flow at Re = 790 Using INS3D and Grid in Fig. 3-4 
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Due to the sparseness of the grid at the symmetry plane boundary, it 
could be argued that the computed structure of the secondary flowfield shown 
in Figs. 3-5 and 3-6 might be grid dependent. To settle this question a 
second calculation at the same Reynolds number was made on the grid in Fig. 
2-2. These results are presented for comparison in Figs. 3-8 through 3-10. 
Figure 3-11 contains results at cross-planes downstream of the turn showing 
the continued development of the swirl pattern. 

The secondary flow pattern for the second case is clearly different 
from that using the first grid. The larger density of grid points near the 
outer wall and symmetry plane intersection has resolved a second counter- 
rotating swirl pattern that appears after mid-turn and grows to be the same 
size as the primary swirl a position near the exit plane of the geometry. 
Because of the resolution of finer details with flow using the grid in Fig. 
2-2, this grid was employed in the turbulent computation for the same 
geometry. 

3 . 2 TURBULENT COMPUTATIONS 


Two-dimensional and axisymmetric turbulent results in the TAD geometry 
using FDNS3D with the standard k-c model are displayed in Figs. 3-12 and 
3-13. No boundary layer separation was observed in either case for a Rey- 
nolds number of one million. Effects of the increase in cross sectional 
area for the axisymmetric case is evident when the two results are contras- 
ted. Each case was run for 1000 iterations on the NASA-MSFC Cray-XMP. 

The geometry in Fig. 2-2 was used to run a turbulent computation at 
Reynolds number of 40,000 using INS3D with the k-e subroutine extracted 
from FDNS3D . This calculation was again run for 1000 iterations. For 
comparison the laminar computation and the turbulent computation Cray-XMP 

-4 

CPU run times were 0.7 x 10 seconds per iteration per node and 2.0 x 
10 4 seconds per iteration per node. Figures 3-14 through 3-20 show 
various aspects of the turbulent square duct computation. The secondary 
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Pressure Contours at Square Duct Mid Plane (Left) and Three-Quarter 
Plane for Laminar Flow at Re = 790 Using INS3D and Grid in Fig. 2-2 


I 


3-10 

LOCKHEED-HUNTSVILLE ENGINEERING CENTER 













LMSC-HEC TR F225988 



I 

I 

I 


Fig. 3-13 


Velocity Vectors (Top) , Velocity Magnitude Contours and Static 
Pressure Contours (Bottom) for Axisymmetric Turbulent Computa- 
tions at Re = 1,000,000 Using FDNS3D with Standard k-e Model 
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Fig. 3-16 Components of Velocity Vectors and Velocity Magnitude Contours in 
Square Duct Cross Planes Corresponding to Two Duct Widths Down- 
stream of Turn (Top) and from Duct Widths Downstream of Turn for 
Turbulent Flow at Re = 40,000 Using INS3D and Standard k-e Model 
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Fig. 3-17 Velocity Vectors and Static Pressure Contours at Symmetry Plane 
(Upper Plots) and Two Planes Above, in a Region Just Downstream 
of 90 deg Turn in Square Duct for Turbulent Computation at 
Re = 40,000 Using INS3D with Standard k-e Model 
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Fig. 3-18 Velocity Vectors and Static Pressure Contours at the Three-Quarter 
Plane (Upper Plots) and Two Planes Below, in a Region Just Down- 
stream of 90 deg Turn in Square Duct for Turbulent Computations at 
Re = 40,000 Using INS3D with Standard k-e Model 
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Pig. 3-19 Components of Velocity Vectors (for Turbulent Computation) 

in Vertical Planes Starting at First Plane Off Inner Wall (Top) 
and Proceeding Eight Planes Toward Mid Channel in Region Just 
Past Mid turn in Square Duct to Three Duct Widths Downstream 
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flow swirl structure evolves in a manner similar to the laminar case with 
the exception that the position of the primary clockwise rotating pattern is 
more tightly fitted into the upper inner wall corner of the duct. This is 
also accompanied by a streamwise adverse pressure gradient on the inner wall 
which is stronger near the symmetry plane than at the upper wall. The 
pressure hill causes the boundary layer to separate for a small region 
downstream of the turn causing a separation bubble on this wall. 

Evidence of this is clearly visible in the plots shown in Figs. 3-17 

through 3-19. In Fig. 3-20, we display the pressure distribution on the 

inner, outer, and upper wall of the duct. 
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4 . CONCLUDING REMARKS 

The Computational Mechanics staff at the Lockheed-Huntsville Engine- 
ering Center have worked closely with NASA-MSFC Computational Fluid Dynamics 
personnel to provide NASA with up-to-date analytical computational results 
and computational tools for SSME engine redesign. Full three-dimensional 
laminar and turbulent SSME/HGM computations have been performed and results 
compared to water flow and air flow measurements (Ref. 1). Experience was 
acquired with two state-of-the-art Incompressible Navier-Stokes codes and 
both have been implemented on the NASA-MSFC supercomputer facility. Results 
of the latter work is reported in this document. Each of the codes, INS3D 
and FDNS3D , yield essentially the same results using identical computational 
grids. The INS3D code is, however, much less storage intensive and is recom- 
mended for large, multi-zone, three-dimensional computations. 
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Appendix A 

INSED k-e TURBULENCE MODEL 
SUBROUTINE LISTING 
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c********************************************************************* 

SUBROUTINE KETURB 

1 ( RENL,VNUT, ITT, ITMAX) 

£***** gg-j. oxMENSIONS ************************************************ 

DIMENSION SS (54, 31, 27, 6) ,VNUT(46656) 

COMMON/XYZQ/Q ( 4 , 46656 ) , X ( 46656 ), Y( 46656 ), Z ( 46656 ) 

COMMON/ FDNS 1/ DK ( 466 56 ) , DE ( 466 56 ) ,TT ( 1 ) , VISE ( 46656) 

COMMON 

1/VAR/UC0( 46656) ,VCO( 46656 ) ,WCO( 46656 ) ,F1< 46656 ) 

1/TRAN/ TJO( 46656) ,DSX( 46656) ,DSY( 46656) ,DSZ( 46656) 

COMMON 

1/PROP/ DEN ( 46656 ) , VISC , DENIN , FLOWIN , DENN1 

1/TUR/ SIGK,SIGE ,CMU , Cl ,C2 ,CMU1 , CMU2 , E ,CK , HINUM , SMNUM , ANV1 ( 7000 ) , 

2 YN( 7000 ) , YN1 ( 7000) ,SINX( 7000) ,SINY(7000) ,SINZ(7000) ,ANW1( 7000) , 

3 IBC ( 7000 ) , JBC( 7000 ) , KBC ( 7000 ) , I ITY( 7000 ) , 

4 GEN (46656) ,MC( 46656) , IJLO( 46656 ) , IITO 
COMMON 

1/COEF/ AP( 46656) , SU ( 466 56 ), SP ( 46656 ), SUK ( 46656 ) , 

2 AE ( 4 66 56 ) , AW (46656) , AN (46656), 

3 AS ( 46656 ), AT (466 56) , AB ( 466 56 ) ,AP0 ( 46656 ) 

COMMON 

1/LIMT/ L , M , LT , MT , LI , L2 , M 1 , M2 , ISWK , ALK , ALVIS , N, N1 , N2 , IG , NT, 

3 ISU, ITU,JSU, JTU , KSU ,-KTU ,CBE,L3,L4,M3,M4,N3,N4, ERRK , ERRDE 

CONVERT INS 3D GEOMETRY TO FDNS3D 


LLMAX=31 
KKMAX* 27 
JJMAX* 54 

NNTOT* LLMAX* KKMAX* JJMAX 

IF ( ITT.EQ. ITMAX) THEN 
IF ( ITMAX. LT. 10 ) GO TO 999 

PRINT *, 'WRITING RESTART QOQ AT ITT = ' , ITT 
DO 692 11=1 , JJMAX 
DO 69 2 J J 2= 1 , LLMAX 
DO 692 KK2=1 , KKMAX 

INOD= I 1+ ( KK2-1 ) * J JMAX+( J J2-1 ) *JJMAX*KKMAX 
SS ( I I , J J2 , KK2 , 1 ) = 1. 

SS ( I I , J J2 , KK2 , 6 ) * 1. 

SS(II,JJ2,KK2,2)» 0(2, INOD) 

SS( II,JJ2,KK2,3)= 0(3, INOD) 

SS(II,JJ2,KK2,4)» 0(4, INOD) 

SS(II,JJ2,KK2,5)= 0(1, INOD )/0.4+0.5* (Q( 2 , INOD) ** 2 
1 + Q ( 3 , INOD ) ** 2+Q ( 4 , INOD) **2 ) 

692 CONTINUE 
FSMACH=0 • 0 
ALPHA=0.0 
TIME* 1 . 

WRITE (23, 190) JJMAX, LLMAX, KKMAX 
WRITE ( 2 3 , 19 5 ) FSMACH , ALPHA , RENL , TIME 

WRITE ( 23, 195) ( ( ( ( SS ( I I , JJ2 , KK2 , N6 ) , 11*1 , JJMAX) , JJ2=1 , LLMAX) , 
1 KK2=1 , KKMAX) , N6=l , 6 ) 

TAGGING K-E & TUR. VISC C 

WRITE( 23, 195) (DK (II), 11=1, NNTOT ) , ( DE( II ) , 11=1 , NNTOT ) 

, (VNUT(II) ,II=l,NNTOT) 
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CLOSE ( 23 ) 

IDB= 1 

IF ( IDB.EQ.l) GO TO 999 

C 

C XYZ FILE C 

C 

PRINT *, 'WRITING RESTART XYZ AT ITT 3 * , ITT 
DO 694 11=1 , JJMAX 
DO 694 JJ2=1 , LLMAX 
DO 694 KK2= 1 , KKMAX 

INOD= I 1+ ( KK2-1 ) *J JMAX+ ( JJ2-1 ) * JJMAX*KKMAX 
SS( II,JJ2,KK2,1)= X(INOD) 

SS(II,JJ2,KK2,2) 3 Y(INOD) 

SS(II,JJ2,KK2,3)= Z(INOD) 

694 CONTINUE 

C OPEN (7, FI LE=' SQD52X.BIN' , STATUS 3 1 NEW' , FORM 3 ’ UNFORMATTED ' ) 

WRITE (2,190) JJMAX, LLMAX, KKMAX 

WRITE ( 2,195) ( ( (SS( I I,JJ2,KK2,1) , 11=1, JJMAX) , JJ2=1 , LLMAX) , 

1 KK2=1, KKMAX) , 

2 ( ( (SS(II,JJ2,KK2,2) ,11=1, JJMAX) ,JJ2=1 , LLMAX) , 

3 KK2=1, KKMAX) , 

4 ( ( <SS(II,JJ2,KK2,3) , I 1=1, JJMAX) ,JJ2*1 , LLMAX) , 

5 KK2=1, KKMAX) 

999 CONTINUE 

PRINT STOP AT KETURB NT=NTMAX & ITT 3 ’, ITT 

STOP 

END IF 

190 FORMAT ( 316 ) 

195 FORMAT (8(E12.5,1X) ) 

C 

DO 197 11=1, JJMAX 
DO 197 JJ2=1, LLMAX 
DO 197 KK2=1, KKMAX 
DO 197 JVAR 3 1 , 4 

I NOD 3 I 1+ ( KK2-1 ) * J JMAX+ ( J J2-1 ) * JJMAX* KKMAX 
SS( II , JJ2,KK2, JVAR) 3 Q(JVAR,INOD) 

197 CONTINUE 

C 

CALL IVA4 ( IDATA , IG , NLIMT , IDUM , 2 , 2 , 1 , 0 ) 

CALL IVA4 ( IMN , JMN , KMN , I SWK ,10,21, 30,6) 

CALL RVA4 ( RDUM , ALK , ALVIS , RDUM , 1000. ,0.2, 0.3,0.) 

c 

C CONSTANTS 

C 

DENN1=0. 1160*22. 876/63. 7/0. 4500 
SIGF 3 0.85 

CALL RVA4 ( DEN IN, VI SC ,SIGU,SIGK, 1.0, 1.0/RENL, 1.0, 1.0) 

CALL RVA4 ( SIGE ,CMU , Cl , C2 , 1.30, 0.09, 1.43, 1.92) 

CALL RVA4 ( E , CK , PI , PI , 9.01069, 0.40000, 3.141592654, 1.000) 

CALL RVA4(SIGK,SIGE,C1,C2, 1.0, 1.15, 1.15, 1.9) 

HINUM-1.E30 
SMNUM*1 . E-30 
CMU1=CMU**0 . 25 
CMU2=CMU**0.75 
C 

c EXAMPLE DATA (90-DEG BEND) 

C 

CALL I VA4 ( L , M , N , LL , JJMAX , LLMAX , KKMAX , 0 ) 

CALL IVA4 (L0,M0,N0,LT, L+l, M+l, N+l , L-l) 

CALL IVA4 ( MT , NT , LI , L2 , M-l, N-l, 0, 0) 
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CALL IVA4(M1,M2,N1,N2, 0, 0, 0, 0) 

CALL IVA4 ( ISU , JSU , KSU , ITU , 2, 2 , 2, LT ) 
CALL IVA4( JTU,KTU,L3,L4, MT , NT, 0, 0) 

CALL IVA4 (M3,M4,N3,N4, 0, 0, 0, 0) 

SQUARE DUCT 90 DEGREE BEND 


WRITE (6,1112) 

1112 FORMAT ( 2 X , * INTO KETURB ##’) 

GET BOUNDARY CONTROL PARAMETERS 

CALL DIRCOS 


INITIALIZE DENSITY/VISCOSITY 


TURINT=0 . 0003 
ALAMDA=0 . 0 1 
INODX= 1 

INODY= 1+(M-1)*L*N 
RLENG=ABS ( Y ( INODX ) - Y ( INODY ) ) 

DO 121 INOD=l , NNTOT 
DEN(INOD) =DENIN 
VNUT ( INOD ) *0 . 0 
VISE ( I NOD) =VISC 

UMEAN2= ( Q ( 2 , INOD ) * * 2+Q ( 3 , INOD) ** 2+Q( 4 , INOD) ** 2 ) 

IF(ITT.GT.l) GO TO 121 

DK( INOD)=0. 5*TURINT*UMEAN2 

DE ( I NOD ) * DK ( INOD ) * * 1 . 5/ ( ALAMDA* RLENG ) 

121 CONTINUE 

SET TURBULENCE QUANTITIES 

IF ( ITT.GT. 1 ) GO TO 123 
DO 122 INOD=l , NNTOT 
IF(MC( INOD) .EQ.0) THEN 

DK ( INOD ) *AMAX1 ( 1 . E-05 , DK ( INOD ) ) 

DE( INOD) =AMAX1( l.E-05,DE( INOD) ) 

ENDIF 

122 CONTINUE 

123 CONTINUE 

CALCULATE GRID TRANSFORMATION COEFFICIENTS 

CALL TRANF ( ITT ) 

C 

DO 800 1=1, L 
DO 800 J«1,M 
DO 800 K=1 , N 

C 

INOD= I+(K-1)*L+(J-1)*L*N 
IP1=MIN0 ( 1+1 , L) 

IM1=MAX0( 1-1,1) 

JP1=MIN0( J+1,M) 

JM1=MAX0 ( J-l , 1 ) 

K?1=MIN0(K+1,N) 

KM1=MAX0(K-1, 1) 

PEW=0 . 5 
PNS = 0 . 5 
PTB=0 . 5 
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C 


c 


800 


IF(I.EQ.l.OR.I.EQ.L) PEW-1.0 
IF(J.EQ.l.OR.J.EQ.M) PNS-1.0 
IF(K.EQ.l.OR.K.EQ.N) PTB-1.0 


I ECC- 
IWCC = 
ICNC = 
ICSC- 
ICCT= 
ICCB- 


I P 1 + ( K- 1 ) 
IM 1 + ( K-l ) 
* (K-l) 


*L+( J-l ) * L*N 
I + ( K-l ) *L+( JP1-1 )*L*N 
I + ( K-l ) *L+( JM1-1 )*L*N 
I +( KP1-1 ) *L+( J-l ) *L*N 
I +(KMl-l)*L+( J-l) *L*N 


P 1-PEW* ( X ( IECC ) -X ( IWCC ) ) 

P2-PNS* ( X( ICNC)-X( ICSC) ) 

P3-PTB* ( X ( ICCT)-X( ICCB ) ) 

Ql-PEW* ( Y( IECC)-Y( IWCC) ) 

Q2=PNS*(Y(ICNC)-Y( ICSC ) ) 

Q3-PTB* ( Y { ICCT) -Y ( ICCB ) ) 

R1=PEW*(Z( IECC) -Z( IWCC) ) 

R2=PNS*(Z(ICNC)-Z(ICSC) ) 
r3=PTB* ( Z ( ICCT ) -Z ( ICCB ) ) 

PTR-1 . 0/ ( P 1 * (Q2*R3-Q3*R2)-P2*(Q1*R3-Q3*R1)+P3*(Q1*R2-Q2*R1 ) ) 
CXC-PTR* (Q2*R3-Q3*R2) 

CYC--PTR* (P2*R3-P3*R2) 

CZC=PTR*(P2*Q3-P3*Q2) 

EXC--PTR* (Q1*R3-Q3*R1 ) 

EYC-PTR* (P1*R3-P3*R1) 

EZC--PTR* ( P1*Q3-P3*Q1 ) 

SXC-PTR* (Ql*R2-Q2*Rl ) 

SYC--PTR* ( P1*R2-P2*R1 ) 

SZC-PTR* ( Pl*Q2-P2*Ql ) 

TJOD-TJO( I NOD ) *DEN ( I NOD ) 

UCO ( INOD ) =T JOD* (Q( 2, INOD ) *CXC+Q ( 3 , INOD ) *C YC+Q ( 4 , INOD ) *C ZC ) 
VCO( INOD ) -TJOD* (0(2, INOD) *EXC+Q( 3 , INOD) *EYC+Q( 4 , INOD ) * EZC ) 
WCO ( INOD) -TJOD* (0(2, INOD) *SXC+Q( 3 , INOD) *SYC+Q( 4 , INOD ) *SZC ) 
CONTINUE 


EVALUATE TURBULENT VISCOSITY 

IF(ITT.GT.l) GO TO 316 
DO 310 INOD-1 , NNTOT 

IF( DK( INOD) .LE. SMNUM ) DK( I NOD) -SMNUM 
I F ( DE ( INOD) .LE. SMNUM) DE ( INOD) -SMNUM 
, IF( DE( INOD) .LE. SMNUM) GO TO 312 
TURVIS=DEN( INOD) *CMU*DK( INOD)** 2/DE ( INOD) +VISC 
GO TO 314 
312 TURVIS-VISC 
314 CONTINUE 

VISE ( INOD) -VISE ( INOD) +ALVIS* ( TURVIS-VISE ( INOD) ) 
310 CONTINUE 
316 CONTINUE 


CALCULATE INLET MASS FLOW RATE 


FLOWIN-O . 0 
AREA- 0.0 
1 = 1 

DO 45 J=2,MT 
DO 45 K-l , N 

INOD- I + ( K-l ) *L+( J-l ) *L*N 
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AREA-AREA+1 .0 
FL0WIN*FL0WIN+UC0( INOD) 

45 CONTINUE 

C WRITE (6, 300) IDATA, AREA , FLOWIN 

ITER S C 

SOLUTION PROCEDURES START 

1 CONTINUE 

*#####*###**#**#*#***#*»*** 

CALL CALCK ( ITT ) 

*#*##*#*#*#*****##***#»*#*# 

«***#**##*#*##****#*»*****# 

CALL CALCE ( ITT ) 

##*##!*##**###**#*#**#####* 

###**#***#*##***#**###**#** 

IF ( ITER.GT.l) GO TO 844 

BOUNDARY CONDITIONS 

DO 15 1=2, L 
DO 15 J=2 ,MT 

INOD1 = I +(J-1) *L*N 

INOD2= I + L+(J-1) *L*N 
0(2, NODI ) =0 ( 2 , INOD2 ) 

0(3, INOD1 ) =Q ( 3 , INOD2 ) 

15 CONTINUE 
C 

r EAST OUT (BASED ON INFLOW MASS FLOW RATE) 

C 

UINC=0 . 0 
FLOW=0 . 0 
DAREA=0 . 0 
DO 50 J=2,MT 
DO 50 K*1,N 

INOD1* LT + ( K-l ) *L+(J-1) *L*N 

INOD2* L +(K-1) *L+(J-1) *L*N 

FLOW=FLOW+UCO( INOD1 ) 

DAREA«DAREA+TJO( INOD2 ) . 67 

50 CONTINUE 

UINC S ( FLOW-FLOWIN) 

DO 51 J=2 ,MT 
DO 51 K=1 , N 

INOD1- LT + ( K-l ) *L+(J-1) *L*N 

INOD2= L + ( K-l ) *L+ ( J-l ) * L*N 

VCO( INOD2 ) =VCO( INOD1 ) 

WCO ( INOD2 ) =WCO( I NODI ) 

51 UCO( INOD2 ) =UCO ( I NODI ) -UINC*TJO( INOD2 ) **0 . 67/DARE A 

C 

DO 52 J=2 , MT 
DO 52 K*2 , NT 


I ECC = 

L 

+(K-1) 

*L+ 

(j-i) 

*L*N 

IWCO 

LT 

+ ( K-l ) 

*L+ 

(j-i) 

*L*N 

ICNC- 

L 

+(K-1) 

*L+ 

(j) 

*L*N 

ICSC = 

L 

+ ( K-l ) 

*L+ 

(J-2) 

*L*N 

ICCT= 

L 

+ (K) 

*L+ 

(J-l) 

*L*N 

ICCB- 

L 

+( K-2 ) 

*L+ 

(J-l) 

*L*N 
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I F ( MC ( I ECC ) . EQ . 0 ) THEN 
P1=1.0*(X(IECC)-X( IWCC) ) 

P2=0 . 5* ( X ( ICNC ) -X ( ICSC ) ) 

P3=0 . 5* ( X( ICCT)-X( ICCB) ) 

01=1. 0*(Y(IECC)-Y( IWCC ) ) 

Q2=0.5*( Y(ICNC)-Y(ICSC) ) 

03=0.5* ( Y( ICCT)-Y( ICCB) ) 

Rl=1.0* ( Z( IECC)-Z( IWCC ) ) 

R2=0.5*(Z(ICNC)-Z( ICSC ) ) 

R3=0. 5* ( Z{ ICCT)-Z( ICCB) ) 

PTR=1.0/(P1*(Q2*R3-Q3*R2)-P2*(Q1*R3-Q3*R1)+P3*(Q1*R2-Q2*R1) ) 
CXC-PTR* (Q2*R3~Q3*R2) 

CYC=-PTR* (P2*R3-P3*R2) 

CZC=PTR* (P2*Q3-P3*Q2) 

EXC=-PTR* (Q1*R3-Q3*R1 ) 

EYC=PTR* (P1*R3-P3*R1 ) 

EZC=-PTR* (Pl*Q3-P3*Ql) 

SXC=PTR* (Q1*R2-Q2*R1 ) 

SYC=-PTR* (P1*R2-P2*R1) 

SZC=PTR* (P1*Q2-P2*Q1 ) 

Q( 2 , IECC ) = UCO( I ECC ) * ( EYC*SZC-SYC*EZC ) + 

1 VCO( IECC) * ( SYC*CZC-CYC*SZC)+ 

2 WCO( IECC) *(CYC*EZC-EYC*CZC) 

Q( 3, IECC)= UCO( IECC ) * ( EZC*SXC-SZC*EXC ) + 

1 VCO(IECC)*(SZC*CXC-CZC*SXC)+ 

2 WCO( IECC ) * ( CZC*EXC-EZC*CXC ) 

Q( 4 , IECC ) = UCO( IECC)*(EXC*SYC-SXC*EYC)+ 

1 VCO( IECC ) * ( SXC*CYC-CXC*SYC ) + 

2 WCO( IECC ) * { CXC*EYC-EXC*CYC ) 

C TT( IECC) =TT( IWCC) 

DK ( I ECC ) *DK ( IWCC ) 

DE( IECC) =DE( IWCC) 

ENDIF 

52 CONTINUE 
844 CONTINUE 

IF ( ITER. GT. 1 ) THEN 
DO 521 J = 2 / MT 
DO 521 K=2 , NT 

IECC= L + ( K-l ) *L+ (J-l) * L*N 
IWCC- LT + ( K-l ) *L+ (J-l) *L*N 
DK( IECC) =DK( IWCC) 

DE ( IECC )=DE( IWCC) 

521 CONTINUE 
END IF 

EVALUATE TURBULENT VISCOSITY 

DO 318 INOD-1 , NNTOT 

IF( DK( INOD) .LE. SMNUM ) DK ( INOD) -SMNUM 
IF(DE(INOD) .LE. SMNUM) DE( INOD) -SMNUM 
IF( DE ( INOD) .LE. SMNUM) GO TO 322 
TURVIS=DEN( INOD) *CMU*DK( INOD)** 2/DE (INOD) +VISC 
GO TO 324 
322 TURVIS=VISC 
324 CONTINUE 

VISE ( INOD) -VISE ( INOD) +ALVIS* ( TURVIS-VISE ( INOD) ) 

VNUT ( INOD ) =TURVIS-VISC 
IF(VNUT( INOD) .LE. 0.0) VNUT ( INOD ) =0 . 0 
318 CONTINUE 

C 
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CONVERGENCE CHECK 


EREXT= 0-01 

NODMN = IMN+ ( KMN-1 ) * L+ ( JMN-1 ) *L*N 

DO 911 1=1, L 

FLOW=0.0 

UMON=0 . 0 

DO 912 J=2,MT 

DO 912 K= 1 , NT 

INOD= I + ( K-l ) *L+( J-l ) *L*N 

FLOW=FLOW+UCO ( INOD) 

UMON=UMON+UCO( INOD) *UCO( INOD) 

912 CONTINUE 
911 CONTINUE 

ERRK=ERRK/( 0. 5*TURINT*FLOW*FLOW) 

IF(ITER.EQ.l) WRITE (6,826) 

826 FORMAT( 10X, ' ITER ' , ' NODMN 4X, ’U’ , 10X, ' DK ’ ,9X, 

. 1 DE' ,9X, ' ERRK ' ,7X, 'VISE' ,7X, ’ VNUT ' , 7X, ’ P ' /) 

WRITE (6,824) ITER , NODMN ,0(2, NODMN ) , DK ( NODMN ) , DE ( NODMN ) 
, ERRK , VISE ( NODMN ) , VNUT ( NODMN ) ,Q( 1 , NODMN ) 

I F ( ITER .GE. 20 .AND. ERRK .GT. 1.E05) GO TO 99 
I F ( ITER .GE. NLIMT .OR. ERRK .LE. EREXT) GO TO 99 
ITER=ITER+1 
ERRK=0 . 0 
GO TO 1 

824 FORMAT ( 2X, ’MONITOR ’ , 2 16 , 7 ( IX, E10 . 4 ) ) 
PRINT OUT SOLUTIONS 


99 CONTINUE 
100 FORMAT ( 1515 ) 

300 FORMAT (1X,I5,6E11.3) 

400 FORMAT (IX) 

500 FORMAT (4I5,3E12.4) 

C****** ***** »»***»»»»»***»***»**** »ti»tt»ii*i»**** ****** ******** 

C 

DO 690 11=1, J JMAX 
DO 690 JJ2=1 , LLMAX 
DO 690 KK2=1 , KKMAX 
DO 690 JVAR= 1,4 

INOD= II+( KK2-1 ) *JJMAX+( JJ2-1 ) *JJMAX*KKMAX 
0( JVAR, INOD) = SS(II,JJ2,KK2, JVAR) 

690 CONTINUE 
RETURN 
END 
C 

cttttttmtttttttttttttttttt 

SUBROUTINE DIRCOS 

CMMMHMttttttttt It ****** 

C 

COMMON/FDNS 1/ DM46656 ) , DE ( 46656 ) ,TT( 1 ) , VISE ( 46656 ) 
COMMON/XYZQ/Q( 4, 46656) ,X( 46656) ,Y( 46656) ,2(46656) 

COMMON 

1/TRAN/ TJO( 46656) ,DSX( 46656 ) ,DSY( 46656 ) ,DSZ< 466 56) 

COMMON 

1/TUR/ SIGK , SIGE , CMU , Cl , C2 , CMU1 , CMU2 , E,CK, HINUM , SMNUM , ANVl ( 7000 ) , 

2 YN ( 7000 ) , YN1 ( 7000 ) ,SINX(7000) ,SINY(7000) ,SINZ(7000) ,ANW1(7000) , 

3 IBC ( 7000 ) , JBC ( 7000 ) , KBC ( 7000 ), IITY( 7000 ) , 

4 GEN (46656) ,MC(46656) , IJLO( 46656 ), IITO 
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1/LIMT/ L,M,LT,MT,L1,L2,M1,M2, ISWK, 

2 ALK,ALVIS,N,N1,N2,IG,NT, 

3 ISU, ITU , JSU , JTU , KSU , KTU ,CBE,L3,L4,M3,M4,N3,N4, 

4 ERRK , ERRDE 


SET DOMAIN BLOCKAGE CONTROL PARAMETER 

-SCALAR BLOCKAGE : MC(INOD)=l 

DO 777 1=1, L 
DO 777 J=1,M 
DO 777 K = 1 , N 

INOD= I + (K-1)*L + ( J-l ) *L*N 
MC ( INOD) =0 

I F ( I . GE . LI . AND .I.LE.L2.AND.J.GE.M1. AND . J . LE . M 2 . AND . 

1 K .GE.N1 . AND.K . LE.N2 ) MC(INOD)=l 

IF( I.GE.L3.AND. I.LE.L4.AND. J . GE . M3 . AND. J . LE .M4 . AND. 

1 K.GE.N3.AND.K.LE.N4) MC(INOD)=l 


FOR SQUARE DUCT 


IF(J.EQ. l.OR. J.EQ.M.OR.K.EQ.N) 
777 CONTINUE 


MC ( INOD ) =1 


CALCULATE BOUNDARY GRID SIZES AND ORIENTATIONS 

1 11 = 1 

DO 30 1=2, LT 
DO 30 J=2 ,MT 
DO 30 K=2 , NT 

INOD= I + ( K-l ) *L + ( J-l ) * L*N 
IF(MC( INOD) .NE. 0) GO TO 30 

— 27 NODE I DENDI FI CATIONS 

CENTER PLANE 


IWCC= 1-1 
I WNC= 1-1 
IWSC= 1-1 
ICCC= I 
ICNC= I 
ICSC= I 
IECC= 1+1 
IENC= 1+1 
IESC= 1+1 


+ ( K-l ) *L +( J-l ) *L*N 
+(K-1)*L + (J) *L*N 
+ ( K-l ) *L +( J-2 ) *L*N 
+ ( K-l ) *L + ( J-l ) *L*N 
+( K-l ) *L +(J) *L*N 
+ ( K-l ) *L + ( J-2 ) *L*N 
+ ( K-l ) *L +( J-l ) *L*N 
+(K-1)*L +(J) *L*N 
+ ( K-l ) *L + ( J-2 ) *L*N 


TOP PLANE 


IWCT= 1-1 + ( K ) 
IWNT= 1-1 + ( K ) 
IWST= 1-1 + ( K ) 
ICCT= I +(K) 

ICNT= I + ( K ) 

ICST= I + ( K ) 

IECT= 1+1 + ( K ) 

IENT= 1+1 + ( K ) 
IEST= 1+1 + ( K ) 


*L + ( J-l ) *L*N 
*L +(J) *L*N 
# L + (J-2 ) *L*N 
*L + ( J-l ) *L*N 
*L +(J) *L*N 
*L + (J-2 ) *L*N 
*L + ( J-l ) *L*N 
*L +(J) *L*N 
*L + ( J-2 ) *L*N 
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C BOTTON PLANE 

IWCB= 1-1 + ( K-2 ) * L + ( J-l ) *L*N 
IWNB= 1-1 + ( K-2 ) *L +(J) *L*N 
IWSB = 1-1 + (K-2 ) *L + ( J-2 ) *L*N 
ICCB= I + ( K-2 ) *L + ( J-l ) *L*N 

ICNB= I + ( K-2 ) * L +(J) *L*N 

ICSB= I + ( K-2 ) * L +( J-2 ) *L*N 
IECB= 1 + 1 + ( K-2 ) * L +( J-l ) *L*N 
IENB= 1 + 1 + ( K-2 ) *L + (J) *L*N 
1ESB= 1+1 + ( K-2 ) *L + ( J-2 ) *L*N 

MCT=MC ( IECC ) +MC ( IWCC)+MC{ ICNC)+MC( ICSC)+ 

1 MC ( I CCT ) +MC ( ICCB ) 

I F ( MCT . EQ. 0) GO TO 30 
IF(MC(ICNC) .EQ. 0) GO TO 2 

-NORTH 

IBC ( I II ) *1 
JBC( III)=J 
KBC ( I I I ) =K 
I ITY ( I I I ) =1 
11=1+1 
12=1-1 
K 1=K+1 
K2=K-1 

IF(MC(IWNC) .EQ.O) 12=1 
IF(MC(IENC) .EQ.O) 11=1 
IF(MC(ICNB) .EQ.O) K2=K 
IF(MC(ICNT) .EQ.O) K1=K 
J1=J+1 
J2=J-1 

INOD3= I + ( K -1) *L+( Jl-1 ) *L*N 

INOD4* I + ( K -1) *L+( J2-1 ) *L*N 

INOD5= I + ( Kl-1 ) *L+( Jl-1 ) *L*N 

INOD6= I + ( K2-1 ) *L+( Jl-1 ) *L*N 

INODC* II + ( K -1) *L+ ( J 1-1 ) # L*N 
INODH* 12 + ( K -1) * L+ (Jl-1 ) *L*N 

P L= ( Y( INOD5 )-Y ( INOD6 ) ) * ( Z ( INODC )-Z ( INODH ) )- 
1 ( Z ( I NOD 5 ) -Z ( INOD6 ) )* ( Y (INODC ) -Y (INODH ) ) 

P2= ( Z ( INOD5 )-Z ( INOD6 ) ) * (X( INODC )-X( INODH ) )- 
1 ( X( INOD5 ) -X( INOD6 ) ) * { 2 ( INODC ) -Z ( INODH ) ) 

P3= ( X( I NOD 5 ) -X ( INOD6 ) )* ( Y ( INODC ) -Y (INODH ) )- 
1 ( Y ( INOD5 ) - Y ( INOD6 ) )*( X (INODC ) -X(INODH ) ) 

PQ=SQRT(P1*P1+P2*P2+P3*P3) 

P4=Pl/PQ 
P5=P2/PO 
P6=P3/PQ 
R1=ABS ( 1 . -P4* * 2 ) 

R2=ABS( l.-P5**2) 

R3=ABS ( 1 . -P6* * 2 ) 

SINX ( III) =SQRT ( R1 ) 

SINY ( I I I ) =SQRT ( R2 ) 

SINZ { I II ) =SQRT ( R3 ) 

Q1=X( INOD) -X( INOD3 ) 

02= Y ( INOD) -Y ( INOD3 ) 

Q3=Z(INOD)-Z( INOD3 ) 
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AA*SQRT( (Q1-P4)**2+(Q2-P5)**2+(Q3-P6)**2) 
CC-1.0 

BB»SQRT(Q1*Q1+Q2*Q2+Q3*Q3) 

COTH* ( BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN{ III )=BB*ABS(COTH) 

Q1 = X( INOD3)-X( IN0D4) 

Q2=Y( I NOD 3 ) - Y( INOD4 ) 

03*2 ( INOD3 ) -Z ( INOD4 ) 

BB=SQRT( 01 *01+02*02+03*03) 

AA=SQRT ( (Q1-P4 ) * * 2+ (Q2-P5 ) ** 2+ (Q3-P6 ) * * 2 ) 
COTH= ( BB*BB+CC*CC-AA*AA)/ ( 2*BB*CC) 

YN1 ( III )=(BB*ABS( COTH) +YN( III) )*0.5 
I JLO ( INOD ) =1 1 1 
III-III+l 
2 CONTINUE 

IF(MC(ICSC) .EQ. 0) GO TO 3 

SOUTH 

IBC ( 1 1 1 ) = 1 
JBC (III) = J 
KBC ( 1 1 1 ) =K 
IITY( III )=2 
11=1+1 
12=1-1 
K1 =K+ 1 
K2=K-1 

I F ( MC ( IWSC ) . EQ . 0 ) 12=1 
I F ( MC ( IESC ) . EQ . 0 ) 11=1 
I F ( MC ( ICSB ) . EQ • 0 ) K2=K 
I F ( MC ( ICST ) . EQ • 0 ) K1=K 
J1=J-1 
J 2= J + 1 

INOD3= I + ( K -1) *L+( Jl-1 ) * L*N 

INOD4= I + ( K -1) *L+(J2-1) *L*N 

INOD5= I + ( Kl— 1 ) *L+ (Jl-1 ) *L*N 

INOD6= I + ( K2-1 ) *L+( Jl-1 ) *L*N 
I NODC= II +(K -1) *L+( Jl— 1 ) *L*N 
INODH= 12 + ( K -1) *L+ ( J 1—1 ) *L*N 

Pl=( Y( INOD5)-Y(INOD6) ) * ( Z ( INODC )-Z ( INODH ) )- 
1 ( Z ( INOD5 ) -Z ( INOD6 ) )*( Y( INODC )-Y( INODH) ) 

P2=(Z(INOD5)-Z(INOD6) )*(X( INODC )-X( INODH) )- 
1 ( X ( I NODS ) -X ( INOD6 ) ) * ( Z ( INODC ) -Z ( INODH ) ) 

P3= ( X( INOD5 )-X( INOD6 ) )*( Y( INODC )-Y( INODH ) )- 
1 ( Y ( INOD5 ) -Y ( INOD6 ) )* (X( INODC)-X( INODH ) ) 

PQ=SQRT(P1*P1+P2*P2+P3*P3) 

P4=P1/PQ 

P5=P2/PQ 

P6=P3/PQ 

R1=ABS(1.-P4**2) 

R2=ABS ( 1 . -P5* * 2 ) 

R3=ABS(1.-P6**2) 

S INX ( I I I ) =SQRT ( Rl ) 

SINY( III) =SQRT ( R2 ) 

SINZ ( I I I ) =SQRT ( R3 ) 

Q1=X ( INOD ) -X ( INOD3 ) 

Q2= Y ( INOD ) -Y ( INOD3 ) 

03=Z( INOD)-Z( INOD3) 


A-10 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 



non 


LMSC-HEC TR F225988 


3 


AA=SQRT( (Q1-P4)**2+(Q2-P5)**2+(Q3-P6)**2) 
CC-1.0 

BB=SQRT(Q1*Q1+Q2*Q2+Q3*Q3) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN ( I I I ) =BB* ABS ( COTH ) 

Q1 =X ( INOD3 ) -X ( IN0D4 ) 

Q2= Y ( I NOD 3 ) - Y( INOD4 ) 

Q3=Z(INOD3)-Z( IN0D4 ) 

BB=SQRT( 01 *01+02*02+03 *03) 

AA=SQRT ( (Q1-P4 ) * * 2+ ( Q2-P5 ) * * 2+ ( Q3-P6 ) **2 ) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN 1 ( 1 1 1 ) = ( BB* ABS ( COTH ) + YN (III) ) * 0 . 5 
IJLO( INOD ) =1 I I 
111=111+1 
CONTINUE 

IF ( MC ( IECC ) .EQ. 0) GO TO 4 


EAST 


C 


C 


IBC ( III )=I 
JBC ( I I I ) =J 
KBC( III )=K 


I ITY ( I I I ) =3 

J1=J+1 

J2=J-1 

K1=K+1 

K2=K-1 

IF(MC(IESC) .EQ.O) 

J2=J 


IF(MC(IENC) .EQ.O) 

J1=J 


IF(MC(IECB) .EQ.O) 

K2=K 


IF(MC(IECT) .EQ.O) 

K1 = K 


11=1+1 

12=1-1 

INODA= 11 +(K1 -1) 

*L+( J-l ) 

*L*N 

INODB= 11 +(K2 -1) 

*L+( J-l) 

*L*N 

INODC= 11 + ( K -1) 

*L+( Jl-1 ) 

*L*N 

INODD= 11 + ( K -1) 

*L+( J2-1 ) 

*L*N 

INODE= 11 + ( K -1) 

*L+(J-1) 

*L*N 

INODF= 12 + ( K -1) 

*L+( J-l ) 

*L*N 


P 1 = ( Y ( INODA ) -Y ( INODB ) ) * ( 2 ( INODC ) -2 ( INODD ) ) - 
1 ( Z ( INODA )-Z ( INODB) )*( Y { INODC ) -Y( INODD) ) 

P2=( Z( INODA )-Z( INODB) )*( X ( INODC ) -X( INODD) )- 
1 (X(INODA)-X( INODB) )*( Z ( INODC ) -2 ( INODD) ) 

P 3= ( X ( INODA ) -X ( INODB ) ) * ( Y ( INODC ) -Y ( INODD ) ) - 
1 ( Y ( INODA) -Y ( INODB ) )* (X( INODC)-X( INODD) ) 

PQ=SQRT(P1*P1+P2*P2+P3*P3) 

P4=P1/PQ 

P5=P2/PQ 

P6=P3/PO 


R1=ABS(1.-P4**2) 
R2=ABS( 1 .-P5**2) 
R3=ABS(1.-P6**2) 
SINX( III ) =SQRT( R1 ) 
SINY ( III) =SQRT ( R2 ) 
SINZ ( III) =SQRT ( R3 ) 
01=X( INOD )-X( INODE) 
Q2=Y ( INOD ) -Y ( INODE) 
Q3=Z(INOD)-Z(INODE) 
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AA=SQRT( (Q1-P4 ) ** 2+ ( Q2-P5 ) ** 2+ ( Q3-P6 ) **2 ) 
CC=1.0 

RB=SQRT(Q1 *01+02*02+03*03) 
COTH-(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN ( 1 1 1 ) =BB*ABS ( COTH ) 

Q1=X( INODE)-X{ INODF ) 

02= Y ( INODE )-Y( INODF) 

Q3= Z ( INODE ) -Z ( INODF ) 
BB=SQRT(Q1*Q1+Q2*Q2+Q3*Q3) 

AA=SQRT( (Q1-P4 ) ** 2+ ( Q2-P5 ) ** 2+ ( Q3-P6 ) ** 2 ) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN1 ( III )=(BB*ABS( COTH) +YN( III) )*0.5 

I JLO ( I NOD ) = I I I 

111=111+1 

CONTINUE 

IF ( MC ( IWCC ) .EQ. 0) GO TO 5 


•WEST 


C 


c 


IBC( III )=I 
JBC ( I I I ) = J 
KBC( III )=K 
I ITY ( I I I ) =4 
J1=J+1 
J2=J- 1 
K1=K+1 
K2 = K- 1 

IF(MC( IWSC) .EQ.O) J 2= J 
IF ( MC ( IWNC) .EQ.O) Jl-J 
IF(MC( IWCB) .EQ.O) K2=K 
IF(MC(IWCT) .EQ.O) K1=K 
11 = 1-1 
12 = 1+1 


INODA= 

11 

+(K1 -1) 

*L+( J-l ) 

*L*N 

INODB= 

11 

+(K2 -1) 

*L+( J-l ) 

*L*N 

INODC = 

11 

+(K -1) 

*L+( Jl-1 ) 

*L*N 

INODD= 

11 

+(K -1) 

*L+ ( J 2-1 ) 

*L*N 

INODE = 

11 

+(K -1) 

*L+( J-l) 

*L*N 

INODF* 

12 

+(K -1) 

*L+( J-l ) 

*L*N 


Pl=( Y{INODA)-Y( INODB) ) 
1 ( Z ( INODA) -Z ( INODB ) ) 

P2= ( Z ( INODA) -Z ( INODB ) ) 
1 ( X( INODA )-X( INODB) ) 

P3= ( X ( INODA) -X ( INODB ) ) 
1 ( Y ( INODA) -Y ( INODB ) ) 

PQ=SQRT(P1*P1+P2*P2+P3 
P4=P1/PQ 
P5=P2/PQ 
P6=P3/PQ 
R1=ABS( l.-P4**2) 


* { Z ( INODC )-Z ( INODD) ) 
*(Y(INODC)-Y( INODD) ) 
* ( X( INODC ) -X( INODD ) ) 
* (Z( INODC )-Z( INODD) ) 
*(Y(INODC)-Y( INODD) ) 
* ( X ( INODC ) -X( INODD) ) 
*P3) 


R2=ABS ( 1 . -P5* *2 ) 
R3=ABS( l.-P6**2) 
SINX( III ) =SQRT( Rl ) 
SINY ( I I I ) “SQRT ( R2 ) 
SINZ ( III ) =SQRT( R3 ) 
Q1=X ( I NOD ) -X ( INODE ) 
Q2= Y( I NOD ) - Y ( INODE ) 
03=Z(IN0D)-Z( INODE ) 
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AA=SQRT( (Q1-P4)**2+(Q2-P5)**2+(Q3-P6)**2) 
CC=1.0 

BB=SQRT(Q1*Q1+Q2*Q2+Q3*Q3 ) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN (1 1 1 ) =BB*ABS ( COTH ) 

Q1=X ( INODE ) -X ( INODF ) 

02= Y ( INODE ) - Y ( INODF ) 

03= Z ( INODE ) -Z ( INODF ) 
BB=SQRT(Q1*Q1+Q2*Q2+Q3*Q3) 

AA=SQRT( (Q1-P4)**2+(Q2-P5)**2+(Q3-P6)**2) 

COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN1( III) =(BB*ABS( COTH )+YN( III) ) *0.5 
IJLO( INOD)=III 
111=111+1 
CONTINUE 

IF ( MC ( ICCT ) .EQ. 0) GO TO 6 


TOP 


C 


C 


IBC( III )=I 
JBC ( 1 1 1 ) = J 
KBC( III )=K 
IITY( 1 1 1 ) = 5 
11 = 1+1 
12 = 1-1 
J1=J+1 
J2=J-1 


IF ( MC ( IWCT ) . EQ . 0 ) 

12=1 


IF(MC(IECT) .EQ.O) 

11 = 1 


IF(MC(ICST) .EQ.O) 

J2=J 


IF ( MC ( ICNT ) .EQ.O) 

J1=J 


K1=K+1 



K2=K-1 



IN0D1= I +(K1 -1) 

*L+( J-l ) 

*L*N 

INOD2= I +(K2 -1) 

*L+( J-l ) 

*L*N 

INOD5= I + ( K 1 — 1 ) 

*L+ ( Jl-1 ) 

*L*N 

INOD7= I + ( Kl-1 ) 

*L+ ( J2-1 ) 

*L*N 

INODA= 11 + { K 1 -1) 

*L+( J-l ) 

*L*N 

INODG= 12 + ( Kl— 1 ) 

*L+ ( J— I ) 

* L*N 


Pl= ( Y ( I NOD 5 ) -Y ( INOD7 ) ) * ( Z ( INODA) -Z ( INODG ) )- 
1 ( Z ( INOD5 ) -Z ( INOD7 ) )*( Y ( INODA) -Y( INODG ) ) 

P2* ( Z ( INOD5 ) -Z ( INOD7 ) ) * ( X( INODA)-X( INODG) )- 
1 ( X( INOD5 )-X( INOD7 ) ) * ( Z ( INODA) -Z ( INODG) ) 

P3= ( X( I NOD 5 ) -X( INOD7 ) )*( Y( INODA)-Y( INODG) )- 
1 ( Y ( INOD5 ) -Y ( INOD7 ) )*( X( INODA ) -X( INODG ) ) 

P0»S0RT(P1*P1+P2*P2+P3*P3) 

P4=P1/PQ 
P5=P2/PQ 
P6=P3/PQ 
R1=ABS (l.-P4**2) 

R2=ABS ( 1 . -P5* * 2 ) 

R3=ABS( l.-P6**2) 

SINX ( I I I ) =SQRT ( R1 ) 

SINY ( I I I ) =SQRT ( R2 ) 

S INZ ( I II ) =SQRT ( R3 ) 


Q1=X( INOD)-X( INODl ) 
02=Y(INOD)-Y(INOD1) 
Q3=Z( INOD)-Z(INODl) 
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AA=SQRT( (Q1-P4)**2+(Q2-P5)**2+(Q3-P6)**2) 

CC=1.0 

BB=SQRT(Ql*Ql+Q2*Q2+Q3*Q3 ) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN ( 1 1 1 ) =BB*ABS ( COTH ) 

Q1=X( IN0D1 ) -X ( IN0D2) 

02= Y ( I NODI ) -Y ( I NOD 2 ) 

03=Z ( I NODI )-Z ( INOD2 ) 

BB=SQRT( 01 *01+02*02+03*03) 

AA=S0RT( (Q1-P4 )**2+(Q2-P5)**2+(Q3-P6)**2) 
COTH=(BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN1( III) = (BB*ABS( COTH )+YN( III) ) *0 . 5 
IJLO( INOD)=III 
111=111+1 
6 CONTINUE 

IF(MC(ICCB) .EO. 0) GO TO 30 
BOTTON 


C 


C 


IBC( III )=I 
JBC ( III )=J 
KBC ( I II ) =K 
IITY( III )=6 
11 = 1+1 
12=1-1 
J1=J+1 
J2=J-1 

IF ( MC ( IWCB ) . EO • 0 ) 1 2=1 
IF ( MC ( IECB ) . EQ . 0 ) 11 = 1 
IF(MC( ICSB) .EQ.O) J2=J 
I F ( MC ( ICNB ) . EQ . 0 ) J1=J 
K 1 =K- 1 
K2=K+1 


IN0D1= I +(K1 -1) 

*L+ ( J-l ) 

*L*N 

INOD2= I +(K2 -1) 

*L+ ( J-l ) 

*L*N 

INOD5= I + ( K 1 — 1 ) 

*L+( Jl-1 ) 

*L*N 

INOD7= I + ( Kl-1 ) 

*L+ ( J2-1 ) 

*L*N 

INODA= 11 +(K1 -1) 

*L+ ( J-l ) 

*L*N 

INODG= 12 + ( K 1— 1 ) 

*L+( J-l ) 

* L*N 


1 

1 

1 


Pl = ( Y ( I NOD 5 ) -Y ( IN0D7 ) ) * ( Z ( INODA) -Z ( INODG ) )- 
( Z ( INOD5 ) -Z ( INOD7 ) ) * ( Y( INODA) -Y( INODG ) ) 
P2= ( Z ( INOD5 )-Z ( INOD7 ) ) * ( X ( INODA ) -X( INODG ) )- 
( X ( INOD5 ) -X( INOD7 ) )*( Z ( INODA) -Z ( INODG) ) 
P3= ( X ( I NOD 5 ) -X ( INOD7 ) ) * ( Y ( INODA)-Y ( INODG) )- 
( Y ( I NODS ) -Y ( INOD7 ) )*( X ( INODA) -X( INODG) ) 
P0=SQRT(P1*P1+P2*P2+P3*P3) 

P4=P1/PQ 

P5=P2/PO 

P6=P3/PO 

R1=ABS(1.-P4**2) 

R2=ABS ( 1 .-P5**2 ) 

R3=ABS(1.-P6**2) 

SINX ( III) =SQRT ( Rl ) 

SINY(III)=S0RT(R2) 

SINZ ( I I I ) =SQRT ( R3 ) 

Q1=X( I NOD ) -X ( INODl ) 

Q2=Y(INOD)-Y(INODl) 

Q3=Z(INOD)-Z(INODl) 
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AA=SQRT( (Q1-P4) # *2+(Q2-P5)**2+(Q3-P6)**2) 

CC«1.0 

BB»SQRT( 01 *01+02*02+03*03) 

C0TH=( BB*BB+CC*CC-AA*AA)/( 2*BB*CC) 

YN ( 1 1 1 ) =BB*ABS ( COTH ) 

Q1=X ( IN0D1 ) -X( IN002 ) 

02= Y ( I NODI ) -Y ( INOD2 ) 

03= Z ( I NODI ) -Z ( INOD2 ) 

BB=SQRT( Q1 *01+02*02+03*03) 

AA=SQRT ( (Q1-P4 ) # *2+(Q2-P5)**2+(Q3-P6)**2) 

C0TH= ( BB*BB+CC*CC-AA*AA)/< 2*BB*CC) 

YN 1 ( I I I ) = ( BB* ABS ( COTH ) + YN ( 1 1 1 ) ) * 0 . 5 
I J LO ( INOD)=III 
111=111+1 
30 CONTINUE 
I ITO= I I I- 1 

C WRITE (6,100) L , M , N , I ITO 

100 FORMAT (415) 

RETURN 

END 

MHHHH ******************* 

SUBROUTINE TRANF ( ITT) 

###**#***##**#*****# 1 *****##* 

COMMON/FDNS1/ DK ( 466 56 ), DE ( 46656 ), TT ( 1 ), VISE ( 466 56 ) 
COMMON/XYZQ/Q( 4 , 46656) ,X( 46656 ) ,Y( 46656 ), Z ( 46656 ) 

COMMON 

l/VAR/UCO( 46656) ,VCO( 46656 ) ,WCO( 46656 ) , FI (46656) 

1/TRAN/ TJO( 46656) ,DSX( 46656 ) ,DSY( 46656) ,DSZ( 46656) 

COMMON 

1/PROP/ DEN (46656) , VISC , DENIN , FLOWIN , DENN1 

1/TUR/ SIGK,SIGE,CMU,C1,C2,CMU1,CMU2,E,CK,HINUM,SMNUM,ANV1(7000) , 

2 YN ( 7000 ) , YN1 ( 7000 ) , SINX ( 7000 ) , SINY ( 7000 ) ,SINZ(7000) ,ANW1(7000) , 

3 IBC( 7000 ) , JBC ( 7000 ) ,KBC( 7000) ,IITY(7000) , 

4 GEN ( 46656 ),MC( 46656 ) ,IJLO( 46656 ) ,IITO 
COMMON 

1/COEF/ AP( 46656) ,SU( 46656 ) ,SP( 46656) ,SUK( 46656) , 

2 AE ( 466 56 ) , AW ( 46656 ) , AN( 46656 ) , 

3 AS ( 46656 ), AT (46656) ,AB( 46656) ,AP0( 46656) 

COMMON 

1/LIMT/ L, M, LT, MT, LI, L2, Ml, M2, ISWK,ALK ,ALVIS , N, N1 , N2 , IG , NT, 

3 ISU,ITU,JSU,JTU,KSU,KTU,CBE,L3,L4,M3,M4,N3,N4,ERRK,ERRDE 

CALCULATE GRID TRANSFORMATION COEFFICIENTS 

DO 40 1=1, L 
DO 40 J=1 ,M 
DO 40 K=1,N 

INOD= I +(K-1)*L +(J-1)*L*N 
IP1=MIN0 ( 1+1 , L) 

IM1=MAX0 (1-1,1) 

JP1=MIN0(J+1,M) 

JM1=MAX0 ( J-l , 1 ) 

KP1=MIN0(K+1,N) 

KM1=MAX0 ( K-l , 1 ) 

PEW=0 . 5 
PNS = 0. 5 
PTB=0 . 5 
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C 


C 


IF ( I .EQ-l.OR.I.EO.L) PEW- 1.0 
IF ( J.EQ. l.OR.J.EO.M) PNS-1.0 
IF(K.EQ.l.OR.K.EQ.N) PTB-1.0 


INOD= 
I ECC = 
IWCC = 
I CNC = 
ICSC- 
ICCT = 
ICCB= 


I + ( K-l ) *L + ( J- 1 ) * L*N 


iriT \ r\ ± f 

IM1+ ( K-l ) 

I + ( K— 1 ) 

I + ( K-l ) 

I + (KP1-1 ) *L+( J-l ) *L*N 
I +(KMl-l)*L+( J-l) *L*N 


*L+( J-l ) * L*N 

*L+( JP1-1 ) *L*N 
*L+( JM1-1 ) *L*N 
*T.+ f .f-1 


P 1 =PEW* ( X ( IECC ) -X ( IWCC ) ) 
P2 = PNS* (X( ICNC)-X( ICSC) ) 
P3=PTB* ( X( ICCT )-X( ICCB) ) 
Q1=PEW* ( Y( IECC) -Y( IWCC) ) 
Q2=PNS*(Y( ICNC ) -Y ( ICSC ) ) 
03 = PTB* ( Y( ICCT )-Y( ICCB) ) 
R1=PEW* ( Z ( IECC ) -Z ( IWCC ) ) 
R2=PNS*(Z( ICNC )-Z( ICSC) ) 
R3=PTB*(Z( ICCT )-Z( ICCB) ) 


PTR«1.0/(P1*(Q2*R3-Q3*R2)-P2*(Q1*R3-Q3*R1)+P3*(Q1*R2-Q2*R1) ) 
AD1=P1 # (Q2*R3-Q3*R2) 

AD2=-P2*(Q1*R3-Q3*R1) 

AD3=P3* (Q1*R2-Q2*R1) 

AD4= ( AD1+AD2+AD3 ) 

IF(AD4.EQ.O.) 

1 WRITE (6,1111) I NOD, Pl,P2,P3/QlfQ2,Q3,Rl,R2 f R3 , ADI , AD2 , AD3 , PTR 
PTR 3 1 . 0/AD4 

T JO { INOD ) =ABS ( 1 . 0/PTR ) 

40 CONTINUE 
1111 FORMAT (1X,3I3,1X,9E10.5/,10X,4E10.5) 

C 


DO 50 1=1, L 
DO 50 J= 1 , M 
DO 50 K= 1 , N 

INOD= I + ( K-l ) *L+ ( J-l ) *L*N 

IF(I.GT.l) THEN 

IWCC= I— 1+(K— 1 ) *L+ ( J— 1 ) *L*N 

DSX( INOD ) =SQRT ( ( X( INOD ) -X( IWCC ) ) ** 2+ ( Y( INOD) -Y ( IWCC ) )**2 
1 + ( Z ( INOD) -Z ( IWCC ) ) **2 ) 

ENDIF 

IF(J.GT.l) THEN 

ICSC= I + ( K-l ) *L+(J-1-1)*L*N 


DSY ( INOD ) =SQRT ( ( X( INOD) -X ( ICSC ) ) **2+ ( Y( INOD) -Y ( ICSC ) )**2 
1 + ( Z ( INOD) -Z ( ICSC ) ) * * 2 ) 

ENDIF 

IF(K.GT.l) THEN 

ICCB= I +(K-1-1)*L+(J-1) *L*N 


DSZ ( INOD) =SQRT( (X( INOD)-X( ICCB) ) * # 2+ ( Y( INOD)-Y ( ICCB ) )**2 
1 + ( Z ( INOD ) -Z ( ICCB ) ) **2 ) 

ENDIF 

50 CONTINUE 
RETURN 
END 

C 

C«#############*#######*#######f#t####tf ##»#»####«#« 
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SUBROUTINE CALCK ( ITT ) 

C*#*|*t»*l*##***#*#*»t»t»*»#ttt»#t*»»»»#*tt#****»»»# 

c 

C0MM0N/FDNS1/ DK ( 46656) , DE ( 46656 ) ,TT ( 1 ), VISE ( 46656 ) 
COMMON/XYZQ/Q( 4,46656) ,X( 46656) ,Y( 46656) ,2(46656) 

COMMON 

1/VAR/UC0( 46656) ,VCO( 466 56) ,WCO( 46656 ), FI ( 46656 ) 

1/TRAN/ TJO( 46656) ,DSX( 46656) , DSY( 46656 ), DSZ ( 46656 ) 

COMMON 

1/PROP/ DEN ( 46656 ) , VISC , DENIN , FLOWIN , DENN1 

1/TUR/ SIGK , SIGE , CMU , Cl ,C2 , CMU1 ,CMU2 , E ,CK, HINUM ,SMNUM , ANV1 (7000 ) , 

2 YN ( 7000 ) , YN1( 7000 ) ,SINX( 7000) , SINY( 7000 ), SINZ ( 7000 ) ,ANW1(7000) , 

3 IBC ( 7000 ) , JBC( 7000 ) ,KBC(7000) ,IITY(7000) , 

4 GEN (46656) ,MC( 46656) , IJLO( 46656 ), IITO 
COMMON 

1/COEF/ AP( 46656) ,SU( 46656) , SP ( 46656 ) ,SUK ( 46656 ) , 

2 AE ( 4 6 6 5 6 ) , AW (46656) , AN (46656) , 

3 AS (46656) , AT ( 46656 ) ,AB( 46656 ) ,AP0( 46656) 

COMMON 

1/LIMT/ L,M,LT,MT,L1 , L2 , Ml , M2 , ISWK , ALK , ALVIS , N , N1 ,N2, IG,NT, 

3 ISU , ITU , JSU , JTU , KSU , KTU , CBE , L3 , L4 , M3 , M4 , N3 , N4 , ERRK , ERRDE 

— TRANSPORT EQUATIONS LINERAI ZATION AND SOLVER 

ERRK*0 . 0 
PI»3. 141592654 

-PRESSURE CORRECTION SOLVER STARTS FROM 10 


CALL 

IVA4 ( IS , IT , JS , JT , 

ISU , 

ITU, 

JSU, 

JTU) 

CALL 

IVA4 ( KS , KT , I I , JX , 

KSU , 

KTU, 

0, 

0) 

CALL 

RVA4 (CE ,CW,CN ,CS , 

0.0, 

0.0, 

0.0, 

0.0) 

CALL 

RVA4 ( CT , CB , PI , P2 , 

0.0, 

0.0, 

0.0, 

0.0) 


U, V, W, TM, K & E EQUATIONS 

DO 22 I=IS-1 , IT+1 
DO 22 J=JS-1 , JT+1 
DO 22 K=KS- 1 , KT+ 1 
INOD = I +(K-1)*L + ( J— 1 ) *L*N 
Fl( INOD)=VISE( I NOD ) /S IGK 
22 CONTINUE 

CALCULATE THE SOURCE TERMS OF TRANSPORT EQNS . 

PI=3. 141592654 
PCC=0.0 

EVALUATE LINK COEFF. AND SOURCE TERMS 

DO 20 I 3 IS , IT 
DO 20 J=JS , JT 
DO 20 K=KS , KT 

— 27 NODE IDENDIFICATIONS 

CENTER PLANE 

IWCC= 1-1 + ( K-l ) *L + ( J-l ) *L*N 
I WNC= 1-1 + ( K-l ) *L +(J) *L*N 
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iwsc= 

INOD= 
ICNC = 
ICSC= 
IECC= 
I ENC= 
IESC = 


1-1 +(K-1)*L 
I + ( K-l ) *L 
I + ( K-l ) *L 
I + ( K-l ) *L 
1+1 + ( K-l ) *L 
1+1 +(K-1)*L 
1+1 +(K-1)*L 


+( J-2)*L*N 
+ ( J-l ) *L*N 
+(J) *L*N 
+( J-2)*L*N 
+( J-l )*L*N 
+(J) *L*N 
+( J-2)*L*N 


TOP PLANE 


I WCT= 1-1 + ( K ) 
IWNT= 1-1 + ( K ) 
IWST= 1-1 + ( K ) 
ICCT= I +(K) 

ICNT= I + ( K ) 

ICST= I + ( K ) 

IECT= 1+1 + ( K ) 

I ENT= 1+1 + ( K ) 
IEST= 1+1 + (K) 


*L + ( J-l ) *L*N 
*L +(J) *L*N 
*L + ( J-2 ) *L*N 
*L + ( J-l ) *L*N 
*L +(J) *L*N 
* L + ( J-2 ) *L*N 
*L + ( J-l ) *L*N 
*L +(J) *L*N 
*L + ( J-2 ) *L*N 


BOTTON PLANE 


IWCB= 1-1 + ( K-2 ) *L + ( J-l ) *L*N 
IWNB= 1-1 + ( K-2 ) *L + (J) *L*N 
IWSB= 1-1 + ( K-2 ) * L + ( J-2 ) *L*N 
ICCB= I + ( K-2 ) *L + ( J-l ) *L*N 
ICNB= I + ( K-2 ) *L +(J) # L*N 
ICSB= I + ( K-2 ) *L + ( J-2 ) *L*N 
IECB= 1+1 + ( K-2 ) *L + ( J-l ) *L*N 
I ENB= 1+1 + ( K-2 ) *L +(J) *L*N 
I ESB= 1+1 + ( K-2 ) *L + ( J-2 ) *L*N 


DENC=1 . 0 

IF ( MC( INOD) . EQ. 1 ) DENC=DENN1 
P1=0 . 5* ( X ( IECC ) -X ( IWCC ) ) 

P2=0.5*(X(ICNC)-X(ICSC) ) 

P3=0 . 5* ( X ( ICCT ) -X( ICCB ) ) 

01=0. 5* (Y( IECC )-Y( IWCC) ) 

Q2=0.5*(Y(ICNC)-Y( I CSC ) ) 

03 = 0 . 5* ( Y( ICCT ) — Y ( ICCB) ) 

R1=0 . 5* ( Z ( IECC )-Z( IWCC) ) 

R2=0.5*(Z(ICNC)-Z(ICSC)) 

R3=0.5*( Z( ICCT) -Z( ICCB) ) 

PTR=1.0/(P1*(Q2*R3-Q3*R2)-P2*(Q1*R3-Q3*R1)+P3*(Q1*R2-02*R1)) 
CXQ=PTR* ( Q2*R3-Q3*R2 ) 

CYQ—PTR* ( P2*R3-P3*R2 ) 

C2Q=PTR*(P2*Q3-P3*Q2) 

EXQ=-PTR*(Q1*R3-Q3*R1) 

EYQ=PTR*(P1*R3-P3*R1) 

EZQ=-PTR* ( P1*Q3-P3*Q1 ) 

SXQ=PTR*(Q1*R2-02*R1) 

SYQ=-PTR*(P1*R2-P2*R1) 

SZQ=PTR* ( P1*Q2-P2*Q1 ) 

GAE=0.5*(F1(IECC)+F1(INOD) ) 

GAW=0 . 5* ( FI ( IWCC ) +F1 ( INOD) ) 

GAN=0 . 5* ( FI ( ICNC ) +F1 ( INOD) ) 

GAS=0 . 5*(F1(ICSC)+F1( INOD) ) 

GAT=0 . 5* ( FI ( ICCT ) +F1 ( INOD) ) 

GAB=0 . 5* ( FI ( ICCB ) +F1 ( INOD) ) 

CE=0. 5* (UC0( INOD) +UC0( IECC) ) *DENC 
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CW=0 . 5* ( UCO ( INOD ) +UCO ( IWCC ) ) *DENC 
CN=0 . 5* ( VCO( INOD) +VCO( ICNC) ) *DENC 
CS=0. 5* (VCO( INOD)+VCO( ICSC) ) *DENC 
CT=0. 5* (WCO( INOD) +WCO( ICCT) )*DENC 
CB=0 . 5* ( WCO( INOD) +WCO( ICCB ) ) *DENC 

FACTEW=0 . 5* ( DSX( INOD ) +DSX ( IECC ) ) 

FACTNS-0 . 5* ( DSY( INOD) +DSY( ICNC) ) 

FACTTB-0. 5*(DSZ( INOD)+DSZ( ICCT) ) 
TXEW-(CXQ**2+CYQ**2+CZQ**2)*TJO(INOD) 
TYNS=(EXQ**2+EYQ**2+EZQ**2) *TJO( INOD) 

TZTB- ( SXQ** 2+SYQ**2+SZQ**2 ) *TJO( INOD) 

TXYN-(CXQ*EXQ+CYQ*EYQ+ 

1 CZQ*EZQ)*0.25*TJO(INOD) 

TYZN= ( EXQ*SXQ+EYQ*SYQ+ 

1 EZQ*SZQ)*0.25*TJO( INOD) 

TZXE- ( SXQ*CXQ+SYQ*CYQ+ 

1 SZQ*CZQ)*0.25*TJO( INOD) 

DDE=GAE*TXEW* FACTEW/DSX ( IECC ) + ( GAN-GAS ) *TXYN+ ( GAT-GAB ) *TZXE 
DDW=GAW*TXEW*FACTEW/DSX( INOD) +( GAS-GAN ) *TXYN+ ( GAB-GAT ) *TZXE 
DDN=GAN*TYNS*FACTNS/DSY ( ICNC ) + ( GAE-GAW) *TXYN+ ( GAT-GAB ) *TYZN 
DDS=GAS*TYNS*FACTNS/DSY ( INOD) + (GAW-GAE ) *TXYN+ ( GAB-GAT ) *TYZN 
DDT=GAT*TZTB*FACTTB/DSZ ( ICCT) +(GAE-GAW) *TZXE+(GAN-GAS ) *TYZN 
DDB=GAB*TZTB* FACTTB/DS Z ( INOD) + ( GAW-GAE ) *TZXE+ ( GAS-GAN ) *TYZN 

-HYBRID SCHEME 

AE ( INOD ) = ( AMAX1 ( ABS ( 0 . 5*CE ) , DDE ) -0 . 5*CE ) 

AW( INOD) =(AMAX1( ABS ( 0.5*CW) , DDW) +0 . 5*CW) 

AN ( I NOD ) = ( AMAX1 ( ABS ( 0 . 5 *CN ) , DDN ) -0 . 5*CN ) 

AS ( INOD) - ( AMAX1 ( ABS ( 0 . 5*CS ) , DDS ) +0 . 5*CS ) 

AT( INOD)- (AMAX1( ABS (0.5*CT) ,DDT)-0. 5*CT) 

AB ( INOD ) * ( AMAX1 ( ABS ( 0 . 5*CB ) , DDB ) +0 . 5*CB ) 


CPG= ( CE-CW+CN-CS+CT-CB ) 

CPO-ABS ( CPG ) 

ANEQ= (GAE+GAN) *TXYN 
ASEQ=- ( GAE+GAS ) *TXYN 
ANWQ=- ( GAW+GAN ) *TXYN 
ASWQ = ( GAW+GAS ) *TXYN 

AETQ- ( GAE+GAT ) *TZXE 
AEBQ=-( GAE+GAB ) *TZXE 
AWTQ=- ( GAW+GAT ) *TZXE 
AWBQ= ( GAW+GAB ) *TZXE 
ANTQ= (GAN+GAT) *TYZN 
ANBQ=- ( GAN+GAB ) *TYZN 
ASTQ=- ( GAS+GAT ) *TYZN 
ASBQ- ( GAS+GAB ) *TYZN 
SU ( INOD) =ANEQ*DK ( IENC ) + 

1 ANWQ*DK ( IWNC ) +ASEQ*DK( IESC ) + 

2 ASWQ*DK ( IWSC ) +AETQ*DK ( IECT ) + 

3 AEBQ*DK( IECB) +AWTQ*DK( IWCT) + 

4 AWBQ*DK( IWCB) + ANTQ * DK ( I CNT ) + 

5 ANBQ*DK( ICNB) +ASTQ*DK ( ICST ) + 

6 ASBQ*DK ( ICSB ) 

SUK ( INOD) =CP0 

UCXI=0.5*(Q(2,IECC)-Q(2,IWCC) ) 
UEDA-0 . 5* ( Q( 2 , ICNC ) -Q ( 2 , ICSC ) ) 
USCI=0.5*(Q(2 f ICCT) -Q( 2 , ICCB) ) 
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VCXI=0.5*(Q(3, IECC )-Q( 3 , IWCC ) ) 
VEDA=0.5*(Q(3,ICNC)-Q(3,ICSC) ) 
VSCI=0.5*(Q(3,ICCT)-Q(3,ICCB) ) 

WCXI=0.5*(Q(4, IECC ) -Q ( 4 , IWCC ) ) 

WEDA=0 . 5* (Q( 4 , ICNC)-Q( 4 , ICSC) ) 

WSCI=0.5*(Q(4,ICCT)-Q(4,ICCB) ) 

UX=UCXI*CXQ+UEDA*EXQ+USCI*SXQ 

UY=UCXI*C YQ+UEDA*EYQ+USCI *SYQ 

UZ=UCXI*CZQ+UEDA*EZQ+USCI*SZQ 

VX=VCXI*CXQ+VEDA*EXQ+VSCI*SXQ 

VY=VCXI‘CYQ+VEDA*EYQ+VSCI*SYQ 

VZ=VCXI*CZQ+VEDA*EZQ+VSCI*SZQ 

WX=WCXI*CXQ+WEDA*EXQ+WSCI*SXQ 

WY=WCXI*C YQ+WEDA* E YQ+WSC I * S YQ 

WZ=WCXI*CZQ+WEDA*EZQ+WSCI*SZQ 

GEN ( I NOD ) =VISE ( I NOD ) * ( ( UY+VX ) “ 2 + ( VZ+WY) ** 2+ ( WX+UZ ) **2+ 
1 2* ( UX*UX+VY*VY+WZ*WZ ) -PCC* 2* (UX+VY+WZ ) **2/3 . ) 

IF (MC ( I NOD) .GT. 0 ) GEN ( INOD) =0 . 0 
20 CONTINUE 


CALCULATE SOURCE TERMS FOR CALCK 

GEN IS CREATED ABOVE DO-LOOP 
K-SOURCE 

DO 55 I“IS , IT 

DO 55 J=JS,JT 

DO 55 K=KS , KT 

INOD- I+(K-1)*L+(J-1)*L*N 

SUTM=GEN( INOD) +0.0 

P1=DEN( I NOD) **2 

SP( INOD)=-SUK( INOD)-TJO( INOD) *CMU*P1*DK( INOD) /VISE ( INOD) 
SU ( INOD)=SUTM*TJO(INOD)+SUK( INOD)*DK( INOD) +SU( INOD) 

55 CONTINUE 

MODIFY WALL BOUNDARY CONDITIONS THRU WALL FUNCTIONS 

I F ( IG .NE. 2) GO TO 410 

EVALUATE WALL BOUNDARY CONDITIONS USING WALL FUNCTIONS 

DO 150 111=1/ IITO 
I=IBC( III) 

J=JBC( III ) 

K=KBC( III) 

K 

INOD= 1+ ( K-l ) *L+ ( J— 1 ) *L*N 

YPLN 1 = DEN (INOD) *SQRT (DK( INOD) )*CMUl*YN( IID/VISC 
IF( YPLNl.LE.il. 65) GO TO 111 

TMULT=DEN ( INOD) ‘SORT ( DK ( INOD) ) *CMUl*CK/ALOG( E*YPLN1 ) 

GO TO 112 

111 TMULT=VISC/YN( III) 

112 TAUN1=-TMULT 

GOTO ( 1, 2,3,4, 5,6) , IITY(III) 

1 CONTINUE 

NORTH 

ICNC= I + (K-l ) *L +(J) *L*N 
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DDU*Q( 2, INOD)-Q( 2, ICNC) 

DDV=Q( 3 , INOD )-Q( 3 , ICNC ) 

DDW=Q(4,INOD)-Q(4,ICNC) 

P1=DDU**2+DDV**2+DDW**2 

GEN( INOD) =P1*TAUN1** 2/VISE ( INOD) 

SU( INOD) =HINUM*SQRT( PI ) *ABS ( TAUN1 ) /DEN ( INOD) /SORT ( 0.09) 
SP( INOD)=-HINUM 
AN ( INOD ) *0 . 0 
GO TO 150 

2 CONTINUE 
C 

C SOUTH 

C 

ICSC= I +(K-1)*L + ( J-2 ) *L*N 
DDU=Q ( 2 , INOD) -Q ( 2 , ICSC ) 

DDV=Q ( 3 , INOD ) -Q ( 3 , ICSC ) 

DDW=Q ( 4 , INOD ) -Q ( 4 , ICSC ) 

P1 = DDU* * 2+DDV** 2+DDW* *2 

GEN ( INOD) =P1*TAUN1** 2/VISE ( INOD) 

SU ( INOD ) =HINUM*SQRT ( PI ) *ABS ( TAUN1 ) /DEN ( INOD) /SORT (0.09) 
SP( INOD)=-HINUM 
AS ( INOD) =0.0 
GO TO 150 

3 CONTINUE 
C 

C EAST 

C 

IECC= 1+1 + ( K-l ) *L + ( J-l ) *L*N 
DDU=Q ( 2 , INOD ) -0 ( 2 , I ECC ) 

DDV=Q ( 3 , INOD ) -Q ( 3 , IECC ) 

DDW=Q ( 4 , INOD) -Q ( 4 , I ECC ) 

P1=DDU**2+DDV**2+DDW**2 

GEN( INOD) =P1*TAUN1** 2/VISE ( INOD) 

SU ( INOD) =HINUM*SQRT (PI ) *ABS (TAUN1 )/DEN( INOD) /SORT (0.09) 
SP( INOD)=-HINUM 
AE ( INOD ) =0 . 0 
GO TO 150 

4 CONTINUE 
C 

C WEST 

C 

IWCC= 1-1 + ( K-l ) *L + ( J-l ) *L*N 
DDU=Q ( 2 , INOD) -Q( 2 , IWCC ) 

DDV=Q ( 3 , INOD ) -Q ( 3 , IWCC ) 

DDW=Q( 4 , INOD)-Q( 4 , IWCC ) 

P1=DDU**2+DDV**2+DDW**2 

GEN( INOD)=Pl*TAUNl**2/VISE( INOD) 

SU ( I NOD ) =HINUM*SQRT ( p 1 ) * ABS ( TAUN 1 ) /DEN ( INOD ) /SORT (0.09) 
SP ( INOD)=-HINUM 
AW( INOD ) “0 . 0 
GO TO 150 

5 CONTINUE 
C 

C TOP 

C 

ICCT= I +(K) *L + ( J-l ) *L*N 
DDU*Q( 2 , INOD)-Q( 2 , ICCT) 

DDV*Q ( 3 , INOD) -0 ( 3 , ICCT ) 

DDW=0 ( 4 , INOD ) -Q ( 4 , I CCT ) 

P1=DDU**2+DDV**2+DDW**2 
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GEN ( I NOD) *P1 *TAUN1* * 2/VISE ( I NOD) 

SU ( I NOD ) = H INUM* SQRT ( P 1 ) * ABS ( TAUN 1 ) /DEN ( INOD ) /SQRT (0.09) 
SP{ INOD)=-HINUM 
AT ( INOD ) =0 . 0 
GO TO 150 
6 CONTINUE 

BOTTOM 

ICCB= I + ( K-2 ) * L + ( J-l ) * L*N 

DDU=Q ( 2 , INOD) -0 ( 2 ,ICCB ) 

DDV=Q ( 3 , INOD ) -Q ( 3 , ICCB ) 

DDW-Q ( 4 , INOD) -Q( 4 , ICCB ) 

P1=DDU* * 2+DDV** 2+DDW** 2 

GEN ( INOD ) = P1 *TAUN1* * 2/VISE ( INOD ) 

SU ( INOD) =HINUM*SQRT( PI ) *ABS ( TAUN1 )/DEN( INOD ) /SQRT ( 0 . 09 ) 
SP( INOD)=-HINUM 
AB( INOD) =0 . 0 
GO TO 150 
150 CONTINUE 

410 CONTINUE 

EAST OUT 

I = IT 

DO 200 J=2,MT 
DO 200 K-2 , NT 

INOD* I + { K-l ) *L + ( J-l ) *L*N 
AE ( INOD) *0 . 0 
200 CONTINUE 

DO 511 I-IS-1 , IT+1 
DO 511 J=JS-1 , JT+1 
DO 511 K=KS-1 , KT+1 
INOD= I + ( K-l ) *L + ( J-l ) *L*N 
511 F 1 ( INOD ) =0 . 0 

LINK COEFF. ASSEMBLY AND BLOCKAGES 

DO 500 I=IS , IT 
DO 500 J=JS,JT 
DO 500 K=KS,KT 
INOD* I + ( K-l ) *L + ( J-l ) *L*N 

ANAB=AE( INOD)+AW( INOD ) +AN ( INOD) +AS ( INOD) +AT( INOD) + 

1 AB( INOD) +AP0( INOD) 

AP ( INOD) *ANAB-SP ( INOD ) 

PDUV=1 . 0 
SCAL= 1 . 0 

I F ( I JLO ( INOD ) . NE . 0 ) SCAL=»SMNUM 
IF ( MC ( INOD) .GE. 1 ) THEN 
AP( INOD)=ALK 
AN ( INOD) -0.0 
AS( INOD) =0.0 
AE ( INOD)=0 . 0 
AW( INOD) =0.0 
AT ( INOD) =0 . 0 
AB( INOD) =0.0 
SU ( INOD ) =DK ( INOD) 

PDUV-0.0 
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ENDIF 

UNDER-RELAXATION 

P1=AMAX1 ( SMNUM , ( AP ( INOD )/0 . 8-ANAB ) ) 
AP ( I NOD ) =AP ( INOD) /ALK 


I ECC = 

I +1 + ( K-l ) 

*L+ ( J-l ) 

*L*N 

IWCC = 

I- 

1+(K-1) 

*L+ ( J-l ) 

*L*N 

ICNC = 

I 

+(K-1) 

*L+ ( J ) 

*L*N 

ICSC-* 

I 

+(K-1) 

*L+( J-2 ) 

*L*N 

ICCT= 

I 

+ (K) 

*L+( J-l ) 

*L*N 

ICCB = 

I 

+ ( K-2 ) 

*L+( J-l ) 

*L*N 


RD=-AP ( I NOD ) * DK ( I NOD ) +AE ( INOD ) * DK ( I ECC ) + 

1 AW( INOD ) *DK ( IWCC ) + AN ( INOD) *DK ( ICNC ) + 

2 AS ( INOD ) *DK ( ICSC ) + AT( INOD) *DK( ICCT)+ AB(INOD)* 

3 DK( ICCB ) +SU ( INOD ) + PDUVM 1 . 0-ALF ) *AP ( INOD) *DK ( INOD ) 

SU( INOD)=RD 

ERRK=ERRK+ABS ( RD*SCAL ) 

500 CONTINUE 

LINEAR EQUATIONS SLOVER 

CALL LINERXd, ISWK , IS , JS , KS , IT , JT , KT , L,M , N , FI ) 

CALCULATE MAXIMUM CORRECTION OF CURRENT ITERATION 

DO 555 I=IS , IT 

DO 555 J=JS,JT 

DO 555 K=KS , KT 

INOD= I + ( K-l ) *L + ( J-l ) *L*N 

DK ( INOD ) a DK ( INOD) +F 1 ( INOD) 

555 CONTINUE 
RETURN 
END 
C 

c#######**##*##*##t*«*##*#*»*###*#»*#**#»*#*##t*####*# 

SUBROUTINE CALCE ( ITT ) 

C#***#««# ####*#***###«#»*##***#*##»*#*##*#»#»»** ****** 

c 

COMMON/FDNS 1/ DK( 46656 ) , DE ( 46656 ) ,TT( 1 ) , VISE ( 46656 ) 
COMMON/XYZQ/Q( 4,46656) , X( 46656 ), Y( 46656 ), Z ( 46656 ) 

COMMON 

1/VAR/UCO ( 46656 ) ,VCO( 46656 ) ,WCO( 46656 ) , FI (46656) 

1/TRAN/ TJO( 46656) ,DSX( 46656 ) ,DSY( 46656 ) ,DSZ( 46656) 

COMMON 

1/PROP/ DEN ( 46656 ) , VISCr DENIN, FLOWIN, DENN1 

1/TUR/ SIGK , SIGE ,CMU , Cl ,C2 , CMU1 ,CMU2 , E,CK , HINUM , SMNUM ,ANV1 ( 7000 ) , 

2 YN( 7000) ,YN1( 7000) , SINX ( 7000 ) ,SINY( 7000 ), SINZ ( 7000 ) ,ANW1(7000) , 

3 IBC ( 7000 ) , JBC ( 7000 ) ,KBC( 7000) ,HTY( 7000) , 

4 GEN( 46656 ) ,MC( 46656 ) ,IJLO( 46656 ),HTO 
COMMON 

1/COEF/ AP( 46656) ,SU( 46656 ) ,SP( 46656 ),SUK( 46656 ) , 

2 AE (46656) , AW (46656) , AN (46656) , 

3 AS ( 46656 ) ( AT( 46656 ) » AB( 46656 ) #AP0( 46656 ) 

COMMON 

1/LIMT/ L,M,LT / MT,L1,L2,M1,M2,ISWK,ALK,ALVIS,N,N1,N2,IG,NT, 

3 ISU, ITU,JSU,JTU,KSU,KTU,CBE,L3,L4,M3,M4,N3,N4,ERRK,ERRDE 

C 
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C TRANSPORT EQUATIONS LINERAI ZATION AND SOLVER 

C 


ERRDE-0.0 

PI-3.141592654 


CALL 

IVA4 ( IS, IT , JS , JT , 

ISU, 

ITU, 

JSU, 

JTU ) 

CALL 

IVA4(KS,KT,II,JX, 

KSU , 

KTU , 

0, 

0) 

CALL 

RVA4(CE,CW,CN,CS, 

0.0, 

0.0, 

0.0, 

0.0) 

CALL 

RVA4 ( CT , CB , PI , P2 , 

0.0, 

0.0, 

0.0, 

0.0) 


E EQUATIONS 

DO 22 I-IS-1 , IT+1 
DO 22 J=JS-1 , JT+1 
DO 22 K=KS- 1 , KT+1 
INOD = I + ( K-l ) *L + ( J-l ) *L*N 
F 1 ( I NOD ) -VISE ( I NOD ) /S IGE 
22 CONTINUE 

EVALUATE LINK COEFF. AND SOURCE TERMS 

DO 20 I=IS , IT 
DO 20 J-JS , JT 
DO 20 K=KS , KT 

— 27 NODE I DENDI FI CATIONS 


CENTER PLANE 



IWCC* 

1-1 

+ ( K-l ) *L 

+ ( J-1)*L*N 

IWNC= 

1-1 

+ ( K-l ) *L 

+(J) *L*N 

IWSC=* 

1-1 

+ ( K-l ) *L 

+ (J-2 ) *L*N 

INOD= 

I 

+ ( K-l ) *L 

+{ J-1)*L*N 

ICNC S 

I 

+ ( K-l ) *L 

+(J) *L*N 

ICSC = 

I 

+ ( K-l ) *L 

+( J-2)*L*N 

I ECC= 

1 + 1 

+(K-1)*L 

+ ( J-l ) *L*N 

I ENC = 

1 + 1 

+ ( K-l ) *L 

+(J) *L*N 

IESC s 

1 + 1 

+ ( K-l ) *L 

+ ( J-2)*L*N 

TOP PLANE 



IWCT= 

1-1 

+ (K) *L 

+ ( J-1)*L*N 

IWNT= 

‘ I- 

1 + ( K ) *] 

L +(J) *L*N 

IWST= 

‘ I- 

1 + ( K ) *1 

L + ( J-2 ) *L*N 

ICCT= 

I 

+ (K) *L 

+ ( J— 1 ) *L*N 

' ICNT= 

I 

+ (K) *L 

+(J) *L*N 

ICST= 

I 

+(K) *L 

+ ( J-2 ) *L*N 

I ECT® 

1 + 1 

+(K) *L 

+( J-1)*L*N 


IENT* 1+1 + ( K ) *L +(J) *L*N 
I EST= 1+1 + ( K ) *L + ( J-2 ) *L*N 

BOTTON PLANE 

IWCB- 1-1 + ( K-2 ) *L + ( J-l ) *L*N 
IWNB- 1-1 + ( K-2 ) *L +(J) *L*N 
IWSB= 1-1 + ( K-2 ) *L + ( J-2 ) *L*N 
ICCB- I + ( K-2 ) *L + ( J-l ) *L*N 
ICNB* I +(K-2)*L +(J) *L*N 
ICSB= I + ( K-2 ) *L + ( J-2 ) *L*N 

IECB- 1+1 + ( K-2 ) *L + ( J-l ) # L*N 
C IENB= 1+1 + ( K-2 ) *L +(J) *L*N 
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I ESB= 1+1 +(K-2)*L + ( J-2 ) *L*N 
DENC-1.0 

IF( MC( I NOD) . EQ. 1 ) DENC-DENN1 
P1=0.5*(X(IECC) -X( IWCC ) ) 

P2=0.5*(X(ICNC)-X( I CSC ) ) 

P3=0 . 5* ( X ( ICCT ) -X ( ICCB ) ) 

Q1=0.5*(Y(IECC)-Y( IWCC ) ) 

Q2=0.5*(Y(ICNC)-Y(ICSC) ) 

03 = 0. 5* (Y( ICCT )-Y( ICCB) ) 

R1=0. 5*( Z( IECC)-Z( IWCC ) ) 

R2 = 0. 5*( Z( ICNC)-Z ( ICSC) ) 

R3=0.5*(Z( ICCT )-Z( ICCB) ) 

PTR=1.0/(P1*(Q2*R3-Q3*R2)-P2*(Q1*R3-Q3*R1)+P3*(Q1*R2-Q2*R1) ) 
CXQ-PTR* (Q2*R3-Q3*R2 ) 

CYQ=-PTR* ( P2*R3-P3*R2 ) 

CZQ-PTR* ( P2*Q3-P3*Q2 ) 

EXQ=-PTR*(Q1*R3-Q3*R1) 

EYQ=PTR*(P1*R3-P3*R1) 

EZQ=-PTR* (P1*Q3-P3*Q1 ) 

SXQ=PTR* (Q1*R2-Q2*R1 ) 

SYQ=-PTR*(P1*R2-P2*R1) 

SZQ-PTR* ( P1*Q2-P2*Q1 ) 

GAE=0 . 5*(F1(IECC)+F1( I NOD ) ) 

GAW=0 . 5 * ( FI ( IWCC ) +F1 ( INOD ) ) 

GAN-0 . 5* ( FI ( ICNC)+F1( I NOD ) ) 

GAS-0. 5* ( FI ( ICSC )+Fl( INOD) ) 

GAT=0. 5* ( FI ( ICCT )+Fl( INOD) ) 

GAB“0 . 5* ( FI ( ICCB ) +F1 ( INOD) ) 

CE=0 . 5* ( UCO( INOD ) +UCO( IECC) )*DENC 
CW=0 . 5* ( UCO ( INOD) +UCO ( IWCC ) ) *DENC 
CN=0 . 5* ( VCO( INOD)+VCO( ICNC) )*DENC 
CS=0 . 5* ( VCO( INOD) +VCO( ICSC ) )*DENC 
CT= 0 . 5 * ( WCO ( I NOD ) +WCO ( ICCT) )*DENC 
CB=0 . 5* ( WCO ( INOD) +WCO( ICCB) ) *DENC 

FACTEW=0 . 5* ( DSX( INOD)+DSX( IECC) ) 

FACTNS=0 . 5* ( DSY{ INOD ) +DSY ( ICNC ) ) 

FACTTB-0 . 5* ( DSZ ( INOD) +DSZ ( ICCT ) ) 

TXEW= ( CXQ* * 2+CYQ* * 2+CZQ* * 2 ) *TJO( INOD) 
TYNS=(EXQ**2+EYQ**2+EZQ**2)*TJO( INOD) 

TZTB= ( SXQ** 2+SYQ** 2+SZQ**2 ) *TJO ( INOD) 

TXYN= ( CXQ*EXQ+CYQ*EYQ+ 

1 CZQ*EZQ)*0.25*TJO( INOD) 

TYZN= ( EXQ*SXQ+EYQ*SYQ+ 

1 EZQ*SZQ)*0.25*TJO( INOD) 

TZXE= ( SXQ*CXQ+SYQ*CYQ+ 

1 SZQ*CZQ)*0.25*TJO( INOD) 

DDE-GAE *TXEW* FACTEW/DSX ( I ECC ) + ( GAN-GAS ) *TXYN+ ( GAT-GAB ) *TZXE 
DDW-GAW*TXEW* FACTEW/DSX ( INOD ) + ( GAS-GAN ) *TXYN+ ( GAB-GAT ) *TZXE 

DDN-GAN*TYNS*FACTNS/DSY ( ICNC ) + (GAE-GAW) *TXYN+( GAT-GAB ) *TYZN 
DDS=GAS*TYNS*FACTNS/DSY ( INOD ) + ( GAW-GAE ) *TXYN+ ( GAB-GAT ) *TYZN 

DDT=GAT*TZTB*FACTTB/DSZ( ICCT) +( GAE-GAW) *TZXE+( GAN-GAS ) *TYZN 

DDB=GAB*TZTB*FACTTB/DSZ ( INOD) + ( GAW-GAE ) *TZXE+ ( GAS-GAN ) *TYZN 

« 

HYBRID SCHEME 

AE( INOD)- (AMAX1(ABS( 0.5*CE) , DDE) -0.5*CE) 

AW ( INOD) » ( AMAX1 ( ABS ( 0 . 5*CW) , DDW)+0 . 5*CW) 

AN ( INOD ) = ( AMAX1 ( ABS ( 0 . 5*CN ) , DDN ) -0 . 5*CN ) 
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AS ( INOD ) * ( AMAXl ( ABS ( 0 . 5*CS ) , DDS ) +0 . 5*CS ) 
AT ( INOD) * ( AMAX1 ( ABS ( 0 . 5*CT ) , DDT)-0 . 5*CT ) 
AB( INOD) 3 (AMAX1(ABS(0.5*CB) , DDB ) +0 . 5*CB ) 


CPG= ( CE-CW+CN-CS+CT-CB ) 

CP0=ABS ( CPG ) 

ANEQ= ( GAE+GAN ) *TXYN 
ASE0=- ( GAE+GAS ) *TXYN 
ANWQ 3 - ( GAW+GAN ) *TXYN 
ASWQ 3 ( GAW+GAS ) *TXYN 
AETQ= (GAE+GAT) *TZXE 
AEBQ 3 - ( GAE+GAB ) *TZXE 
AWTQ— ( GAW+GAT ) *TZXE 
AWBQ 3 ( GAW+GAB ) *TZXE 
ANTQ 3 ( GAN+GAT ) *TYZN 
ANBQ=- ( GAN+GAB ) *TYZN 
ASTQ 3 - ( GAS+GAT) *TYZN 
ASBQ= ( GAS+GAB ) *TYZN 
SU ( INOD ) =ANEQ*DE ( I ENC ) + 

1 ANWQ*DE( IWNC)+ASEO*DE( IESC)+ 

2 ASWQ*DE ( IWSC ) +AETQ*DE ( IECT ) + 

3 AEBQ*DE ( IECB ) +AWTQ # DE ( IWCT ) + 

4 AWBQ*DE( IWCB)+ANTQ*DE( ICNT)+ 

5 ANBQ*DE( ICNB)+ASTO*DE(ICST)+ 

6 ASBQ*DE ( ICSB ) 

SUK ( INOD ) -CPO 

20 CONTINUE 


CALCULATE SOURCE TERMS FOR CALCE 


PI=3. 141592654 
PCC-0.0 

E-SOURCE 


DO 65 I 3 IS , IT 
DO 65 J 3 JS, JT 
DO 65 K=KS,KT 

INOD 3 I + ( K-l ) *L + ( J— 1 ) *L*N 
TMDE=DE ( INOD) 

IF ( TMDE . LE . 1 . E-30 ) TMDE-l.E-30 
SUTM 3 C1*CMU*GEN( INOD) *P1*DK( INOD)/ 

1 VISE ( INOD) +0 . 0 

TMDK=DK ( INOD) +SMNUM 

SP( INOD) —SUK ( INOD) -TJO( INOD) *C2*DEN( INOD) *DE( INOD) /TMDK 
SU( INOD) -SUTM*TJO( INOD) +SUK{ INOD) *DE( INOD) +SU( INOD) 

65 CONTINUE 


■MODIFY WALL BOUNDARY CONDITIONS THRU WALL FUNCTIONS 
I F ( IG .NE. 2) GO TO 410 

■EVALUATE WALL BOUNDARY CONDITIONS USING WALL FUNCTIONS 

DO 150 III 3 l,IITO 
I 3 IBC(III) 

J=JBC( III) 

K=KBC (III) 
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INOD* I + ( K-l ) * L + ( J-l ) *L*N 

NORTH/SOUTH ETC ARE SAME FOR E 

E 

TERM=CMU2/ ( CK* YN ( 1 1 1 ) ) 

SU ( I NOD) *HINUM*TERM*DK ( I NOD ) **1 . 5 
SP( INOD)=-HINUM 
150 CONTINUE 

410 CONTINUE 

EAST OUT 

I*IT 

DO 200 J=2,MT 
DO 200 K*2 , NT 

INOD* I + ( K-l ) *L + ( J-l ) *L*N 
AE ( INOD) =0 . 0 
200 CONTINUE 

DO 501 I*IS-1,IT+1 
DO 501 J=JS-1,JT+1 
DO 501 K=KS-1 , KT+1 
INOD* I + ( K-l ) *L + ( J-l ) *L*N 
501 FI ( INOD) =0 . 0 

LINK COEFF. ASSEMBLY AND BLOCKAGES 


DO 500 I=IS,IT 
DO 500 J-JS,JT 
DO 500 K*KS,KT 


INOD* 

I 

+ ( K-l ) 

*L + ( J-l ) 

*L*N 

I ECC= 

1+1+ ( K-l ) 

*L+( J-l ) 

*L*N 

IWCC= 

I- 

1 + ( K-l ) 

*L+( J-l) 

*L*N 

ICNC* 

I 

+ ( K-l ) 

*L+ ( J ) 

*L*N 

ICSC* 

I 

+ ( K-l ) 

*L+( J-2) 

*L*N 

ICCT= 

I 

+ ( K ) 

*L+ ( J-l ) 

*L*N 

ICCB* 

I 

+(K-2) 

*L+ ( J-l ) 

*L*N 


AN AB= AE ( I NOD ) +AW ( I NOD ) +AN (INOD) +AS ( I NOD ) +AT ( INOD ) + 
. AB ( INOD )+AP0( INOD) 

AP ( INOD ) “ANAB-SP ( INOD ) 

PDUV-1.0 
SCAL=1 . 0 

IF( IJLO( INOD) .NE.0) SCAL=SMNUM 
IF ( MC ( INOD) .GE . 1 ) THEN 
AP( INOD) a ALK 
AN ( INOD) *0 . 0 
AS( INOD) *0.0 
AE ( INOD) *0 . 0 
AW( INOD)“0. 0 
AT( INOD)=0.0 
AB( INOD) =0.0 
SU ( INOD) *DE ( INOD) 

PDUV=0 . 0 
ENDIF 

-UNDER-RELAXATION 
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P1=AMAX1 (SMNUM, (AP( INOD)/0 . 8-ANAB) ) 

AP( INOD) a AP( INOD)/ALK 

RD=-AP ( I NOD) *DE ( I NOD) +AE ( I NOD) *DE ( IECC ) +AW( INOD ) 

1 *DE ( IWCC ) +AN ( INOD) *DE ( ICNC ) +AS ( INOD) *DE( ICSC ) 

2 +AT( INOD) *DE ( ICCT )+AB( INOD)*DE( ICCB) 

3 +PDUV* ( 1 4 0-ALDE ) *AP ( INOD) *DE ( INOD) +SU ( INOD) 

SU( INOD)=RD 

ERRDE=ERRDE+ABS ( RD*SCAL ) 

C 

500 CONTINUE 

C 

c LINEAR EQUATIONS SLOVER 

C 

CALL LINERXd , ISWK, IS, JS,KS,IT, JT, KT, L,M ,N ,F1 ) 

C 

C CALCULATE MAXIMUM CORRECTION OF CURRENT ITERATION 

C 

DO 555 I=IS , IT 
DO 555 J=JS,JT 
DO 555 K*KS,KT 

INOD= I + ( K-l ) * L +(J-1) *L*N 

DE ( INOD) =DE ( INOD) +F1 ( INOD) 

555 CONTINUE 
RETURN 
END 
C 

c«# «**#*##«#****»**#*»##***#*###»*******«»******#»#**##* 

SUBROUTINE IVA4 ( IA, IB, IC, ID, JA,JB,JC,JD) 

C### **#****««**# *##*#* t#»»****#t##t ****»»»*»*#*##**«#*** 

c 

IA= JA 
IB= JB 
IC=JC 
ID=JD 
RETURN 
END 
C 

c*#***m*#**#«*i*#*****»#*#**«###t»##*m*t*»t ##«*#§##« 

SUBROUTINE RVA4 ( PA, PB, PC, PD, QA,QB,QC,QD) 

C**##**#****#«t*M***#*t»**###***t****»i**t»*»»«**»»***t 

c 

PA=QA 

PB=QB 

PC=QC 

PD“QD 

RETURN 

END 

C 

C# ########## t #*#*#*****»*#«**«*#******#»#***#**#*** ##*»*#*##»» 

SUBROUTINE LINERX ( ISOL , ISWF , IS , JS, KS , IT, JT,KT, L,M ,N, F ) 

C###*#*«i##**#i#****************#*******#**#**#*****#*t»***#*« 

c 

DIMENSION A( 90 ) , B( 90 ) ,C( 90 ) , D( 90 ) , F( 46656 ) 

COMMON 

1/COEF/ AP( 46656) ,SU( 46656 ) ,SP( 46656 ) ,SUK( 46656 ) , 

2 AE (46656) , AW (46656) , AN (46656) , 

3 AS (46656) , AT ( 46656 ) ,AB( 46656 ) ,AP0( 46656) 

C 

C LINE-RELAXATION USING TDMA 

C 
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DO 90 ISW=1 , I SWF 
ERRF=0.0 

C 

A( JS-1 ) =0 . 0 
DO 120 I=IS,IT 
DO 120 K=KS,KT 

INODJS= I + ( K-l ) *L + ( JS-2 ) *L*N 
C( JS-1 )=F( INODJS) 

DO 121 J=JS,JT 

C 


INOD= 

I 

+(K-1) 

*L +(J-1) 

*L*N 

I ECO 

i+l + Oc-l ) 

*L+(J-1) 

*L*N 

IWCC* 

I- 

1+ ( K-l ) 

*L+( J-l ) 

*L*N 

ICCT* 

I 

+ (K) 

*L+(J-1) 

*L*N 

ICCB* 

I 

+(K-2) 

*L+(J-1) 

*L*N 


A( J ) “AN ( I NOD ) 

B ( J ) *AS ( INOD ) 

C( J)»SU( INOD )+AE( INOD) *F(IECC)+AW( INOD )*F(IWCC)+ 
1 AT ( INOD ) *F ( ICCT ) +AB ( INOD) *F ( ICCB ) 

D( J)=AP( INOD) 

TERM=1.0/(D(J)-B(J)*A(J-1) ) 

A ( J ) =A( J ) *TERM 

121 C(J)»(C(J)+B(J)*C(J-1) ) *TERM 
DO 122 JX=JS,JT 

J= JS+JT—JX 

INOD= I + ( K-l ) *L +(J-1) *L*N 
ICNC= I + ( K-l ) *L+ (J) *L*N 

122 F ( INOD ) =A( J ) * F ( ICNC ) +C ( J ) 

120 CONTINUE 

C 

A( KS-1 ) *0 . 0 
DO 110 I=IS , IT 
DO 110 J=JS,JT 

INODKS= I + ( KS-2 ) *L +(J-1) *L*N 
C ( KS-1 ) = F ( INODKS ) 

DO 111 K=KS , XT 

C 


I NOD = 

I 

+ ( K— 1 ) 

* L +(J-1) 

*L*N 

I ECC* 

1+1+ ( K-l ) 

*L+( J-l) 

*L*N 

IWCC* 

I- 

1+(K-1) 

*L+( J-l ) 

*L*N 

ICNC* 

I 

+(K-1) 

*L+ ( J ) 

*L*N 

ICSC* 

I 

+(K-1) 

*L+( J— 2) 

*L*N 


A(K)=AT( INOD) 

B ( K ) =AB ( INOD) 

C ( K ) =SU ( INOD) +AE ( INOD) *F( IECC)+AW( INOD ) *F ( IWCC ) + 
1 AN ( INOD ) * F ( ICNC ) +AS ( INOD ) * F ( ICSC ) 

D(K)=AP( INOD) 

TERM=1 . 0/ ( D( K ) -B (K)*A(K— 1) ) 

A(K)=A(K)*TERM 

111 C(K)=(C(K)+B(K)*C(K-1) )*TERM 
DO 112 KK-KS,KT 
K=KS+KT-KK 

INOD- I + ( K-l ) *L +(J-1) *L*N 
ICCT= I + ( K ) *L+(J-1) *L*N 

112 F ( INOD ) =A( K ) *F ( ICCT ) +C ( K ) 

110 CONTINUE 

C 

A( IS-1 ) = 0. 0 
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DO 130 J=JS , JT 
DO 130 K=KS , KT 

INODIS 3 IS-1 + ( K-l ) *L +(J-1) * L*N 
C ( IS-1 ) a F ( I NODI S ) 

DO 131 I 3 IS , IT 

C 


INOD= 

I 

+( K-l > 

* L +(J-1) 

*L*N 

IECC= 

1+1+ ( K-l ) 

*L+( J-l ) 

*L*N 

IWCC = 

I- 

1+ ( K-l ) 

*L+ ( J-l ) 

*L*N 

ICNC* 5 

I 

+ ( K-l ) 

*L+(J) 

*L*N 

ICSC = 

I 

+( K-l ) 

*L+( J-2) 

*L*N 

ICCT= 

I 

+ (K) 

*L+ ( J-l ) 

*L*N 

ICCB= 

I 

+ ( K — 2 ) 

*L+( J-l ) 

*L*N 


A ( I ) =AE ( INOD ) 

B ( I ) =AW ( INOD ) 

C ( I ) 3 SU ( INOD ) +AN ( INOD ) *F ( ICNC ) +AS ( INOD) *F( ICSC)+ 

1 AT ( INOD ) * F ( I CCT ) +AB ( I NOD ) * F ( ICCB ) 

D ( I ) *AP ( INOD ) 

ERRF-ERRF+ABS (C(I)-D(I)*F( INOD ) +A( I ) *F ( IECC )+B( I ) *F( IWCC ) ) 
TERM=1.0/(D(I)-B(I)*A(I-1) ) 

A( I ) =A( I ) *TERM 

131 C(I) 3 (C(I)+B(I)*C(I-1) )*TERM 
DO 132 II=IS , IT 

I=IS+IT-I I 

INOD= I + ( K-l ) *L +(J-1) *L*N 
IECC= I+1+(K-1) *L+ ( J-l ) *L*N 

F ( INOD) =A( I ) *F ( IECC ) +C ( I ) 

132 CONTINUE 
130 CONTINUE 

IF ( ISW.EQ.l) ERRF1*AMAX1(ERRF, 0.000001) 

IF ( (ERRF/ERRF1) .LE.l.E-01) RETURN 
90 CONTINUE 
RETURN 
END 
C 
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